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1.0  Summary 

The  principle  goal  of  S  AIC's  Highly  Controlled  Dielectric  Emissivity  (HIDE)  program 
was  to  apply  the  techniques  of  large  area,  low  cost,  continuous  web  manufacturing  to 
fabricate  photonic  band  gap  structures  that  exhibit  finite  width  pass-  and  stop-band 
properties  in  the  8-12  pm  wavelength  regime.  The  3  major  thrusts  of  the  program  were: 

(1)  Develop  web  compatible  techniques  for  manufacturing  single  layer  periodic 
structures 

(2)  Develop  techniques  that  can  be  used  to  stack,  or  layer,  the  single  layer  structures 

(3)  Develop  a  theoretical  framework  for  understanding  the  performance  of  the 
fabricated  structures 

We  developed  several  innovative  techniques  to  fabricate  single  layer  structures.  All  of 
these  approaches  centered  on  the  use  of  microembossed  structures  which  can  be 
fabricated  over  very  large  areas  at  low  cost.  The  most  important  technique,  direct 
transfer,  can  be  used  to  fabricate  3  dimensional  structures  (we  used  hemispheres)  in  one 
substrate  and  transfer  them  to  another  surface. 

We  also  developed  techniques  for  layering  structures.  We  were  able  to  fabricate 
structures  of  up  to  10  layers.  We  concluded  that  rotational  alignment  is  possible, 
however,  we  were  not  able  to  successfully  achieve  translational  alignment. 

Finally,  we  developed  a  semi-analytic  integral  equation  approach  for  the  computation  of 
the  reflection  and  transmission  coefficients  of  a  periodic  single  or  multilayer  array  of 
dielectric  or  metallic  spheres  inside  a  host  medium,  which  may  be  lossy. 

The  program  was  divided  into  3  phases.  Phase  1  (10/97  -  5/98)  has  been  completed  and 
focused  on  developing  the  appropriate  microembossing,  and  liftoff  procedures  to 
fabricate  a  single  layer  with  the  proper  periodicity.  Phase  2  (5/98  -  5/99)  focused  on 
developing  tools  and  techniques  multilayer  fabrication.  In  Phase  3  (6/98  -  10/99),  we 
continued  development  on  multilayer  techniques  to  produce  multilayer  photonic  crystal 
structures.  The  report  below  summarizes  our  fabrication  and  modeling  efforts. 

The  participating  scientific  personnel  is  listed  below: 

SAIC 

Dr.  Albert  Green 
Dr.  Alan  Phillips 
Dr.  Jill  Inman 
Dr.  Convers  Wyeth 
Dr.  Robb  White 


Dr.  Dimos  Dialetes 
Dr.  David  Chemin 
Mr.  Jon  Ward 
Mr.  Steve  Krafsig 

Polaroid/MicroContinuum/Strategic  Data  Resources 

Dr.  Steve  Hawley 
Dr.  Dennis  Slafer 

Howard  University 

Dr.  Clayton  Bates 

NIST 

Dr.  Erich  Grossman 
Dr.  Carl  Reinsemer 

2.0  Fabrication 

Our  basic  approach  involves  single  layer  fabrication  of  periodic  dielectric  structures  that 
can  be  layered  to  form  a  photonic  crystal.  The  ultimate  goal  is  to  produce  materials  using 
our  continuous  web  techniques  that  exhibit  a  photonic  bandgap  in  the  8-12  pm  (1250  - 
800  cm'1)  regime. 

Preliminary  work  revealed  a  somewhat  difficult  materials  issue,  namely  most  of  the 
conventional  materials  used  in  embossing  technology  such  as  polycarbonate  and 
CAB  (cellulose  acetate  butyrate)  have  significant  absorbing  properties  in  the  8  - 12  pm 
region.  After  an  extensive  literature  and  database  search,  we  identified  two  materials  as 
suitable  matrix  materials  for  the  final  structures,  polyethylene  and  a  material  known  as 
Poly  IR2  which  is  a  proprietary  polymeric  product  of  Fresnel  Technologies,  Inc  and  is 
used  to  fabricate  Fresnel  Lenses  for  8  - 12  pm.  These  materials  are  thermally  embossable 
and  moldable.  Unlike  polycarbonate  and  CAB,  solvent  embossing  with  conventional 
solvents  at  room  temperature  did  not  seem  possible  because  of  solubility  limits.  It  should 
be  pointed  out  that  the  availability  of  soluble  polymer  is  not  essential  for  achieving  the 
ultimate  goals  of  the  program.  For  purposes  of  developing  evaluation  prototypes, 
however,  the  ability  to  use  cast  films  that  are  spun  or  rod  coated  would  provide  an  extra 
degree  of  freedom  and  considerably  more  flexibility  in  the  planarization  step. 

Early  in  the  program  we  evaluated  the  materials  issue  more  carefully  as  part  of 
developing  the  fabrication  process.  The  principal  purpose  of  this  evaluation  was  to 
determine  the  range  of  polymeric  materials  that  do  not  have  significant  absorbance  in  the 
wavelength  region  of  interest  and,  based  on  previous  experience,  can  be  fabricated  in  a 
web  compatible  process.  These  considerations  drove  many  of  the  fabrication  approach 
decisions  throughout  the  program  and  are  worth  a  separate  discussion  below. 


2.1  Substrate  Infrared  Spectra 


For  HIDE  fabrication,  a  principle  considerations  were  infrared  spectra  vis  a  vis  the 
absorption  in  8-12  pm  and  the  solubility.  The  results  are  summarized  below. 

There  are  a  number  of  commercially  available  libraries  either  in  printed  or  electronic 
form  that  provides  infrared  absorption  characteristics  of  materials.  The  purpose  of  these 
databases  is  primarily  for  compound  identification  via  the  various  absorption  peaks.  As 
such,  the  listed  spectra  are  normalized  to  absorption  peak  height  and  not  material 
thickness  so  do  not  provide  an  absolute  measure  of  optical  absorbance.  Nonetheless, 
since  we  were  interested  in  materials  which  have  no  peaks  in  the  8-12  pm  range,  we 
believed  the  database  is  a  useful  guide. 

Unfortunately,  the  cost  of  ownership  of  such  a  database  was  much  higher  than  is  justified 
by  this  investigation.  It  was  possible;  however,  to  convince  Nicolet  Instruments  who  are 
purveyors  of  the  electronic  form  of  these  databases  to  conduct  a  one  time,  no-charge 
search  of  the  databases  most  germane  to  our  goals.  The  search  criterion  was  simply  no 
absorption  peaks  in  8  - 12  pm  region.  Three  databases  were  searched  representing  about 
15,000  organic  compounds.  A  complete  listing  of  the  compounds  searched  will  be 
submitted  in  the  final  report.  12  compounds  were  listed  as  the  best  fit  to  the  search 
criteria.  The  following  compounds  were  considered  the  best  candidates  for  web 
compatible  fabrication  consistent  with  the  HIDE  goals. 

Poly(ethylene)  r\  =  450  cps(viscosity) 

Poly(ethylene)  r|  =  6000  cps 
Poly(l-hexadecene)  isotatic 

Poly(ethylene :propylene : diene)  70%  ethylene,  4%  diene 
Poly(styrene:butadiene)ABA  block  (28%  styrene) 

Poly(ethylene:acrylic  acid)  10%  acrylic  acid 
Poly(ethylene:propylene:diene)  70%  ethylene,  4%  diene(repeated) 

Poly(ethylene:acrylic  acid)  10%  acrylic  acid(repeated) 

Paraffin  wax  #1 

Poly(l  -butene)  isotactic,  high  MW 
Versamine  551 
Thixtrol  SR- 100 

Of  these  compounds,  the  last  two,  versamine  and  thixatrol  are  polymer  additives. 
Versamine  is  a  product  of  the  Henkel  Corporation  and  is  a  curing  agent  for  epoxy 
coatings.  Thixtrol  is  a  product  of  Rheox,  Inc  and  is  a  thixatropic  additive  to  control 
viscosity  of  liquid  materials  like  paint. 

With  the  exception  of  the  styrene:butadiene  copolymer,  the  remainder  all  fall  into  the 
general  family  called  polyolefins  which  include  polyethylene  and  polypropylene.  The 
distinguishing  feature  of  these  compounds  with  regard  to  their  IR  spectra  is  that 
regardless  of  whether  the  polymer  is  linear  or  branched,  the  entire  molecule  contains  only 
two  chemical  bonds,  C-C  or  C-H.  The  result  is  absorption  bands  at  3.5, 7,  and  14p. 
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Paraffin  which  is  a  short  chain  polyolefin  with  a  chain  length  of  about  15-40  carbons 
shows  the  typical  IR  absorption  spectra  for  this  kind  of  material.  An  excellent 
introduction  to  the  properties  of  polyolefins  can  be  found  at 

http://www.umr.edu/~wlf/CHEM381/chap22.html 

We  measured  the  transmission  of  a  number  of  substrates  to  determine  their  suitability  for 
use  in  the  program.  Figure  2.1.1  shows  the  transmission  spectra  of  2.0  mils  of 
polyethylene. 

2.2  Substrate  Solubility 

We  focused  on  Polyethylene  due  to  its  low  absorption  in  8-  12pm.  Polyethylene  is 
relatively  insoluble  in  most  solvents,  in  fact  it  is  sometimes  used  to  make  gasoline 
containers.  The  literature  claims  that  it  has  some  solubility  in  hot  o-xylene  and  CS2.  We 
performed  experiments  that  showed  very  little  indication  of  solubility  with  a  variety  of 
solvents  at  room  temperature.  A  possible  alternative  is  paraffin.  We  have  obtained  the 
highest  melting  grade  from  Aldrich.  This  material  is  a  hard,  brittle  material  with  a  milky 
appearance.  We  were  able  to  produce  saturated  solutions  by  dissolving  this  material  in 
petroleum  distillates  such  as  petroleum  ether  at  60°C.  At  this  point  we  believe  that 
normal  polyethylene  can  be  softened  with  the  same  solvents  at  elevated  temperatures. 

From  this  analysis,  we  concluded  that  in  general,  there  is  a  very  limited  selection  of 
materials  that  can  be  used  for  the  lattice  element  in  these  structures  and  in  particular 
organic  materials.  Simple  polyethylene  film  has  the  correct  absorbing  properties  but 
suffers  from  limiting  mechanical  properties.  There  are  basically  two  kinds  of 
polyethylene  high-density  polyethylene  (HDPE)  and  the  low-density  variant  (LDPE).  The 
difference  is  the  result  of  the  degree  of  polymer  branching  cf.  Appendix  I.  LDPE  is 
relatively  soft.  We  found  that  it  was  very  difficult  to  main  reasonable  dimensional 
integrity  following  the  embossing  process.  Low  density  is  more  rigid  but  also  said  to 
suffer  from  shrinkage  associated  with  molding  although  we  have  not  tried  to  emboss  this 
material. 

Although  we  do  not  know  the  exact  composition  of  Poly  IR2  it  is  reasonable  to  assume  it 
is  also  a  variant  of  polyethylene.  We  obtained  Poly  IR2  samples  and  found  them  to  be 
waxy  and  milky  appearance.  We  presumed  that  is  probably  a  polyethylene  copolymer 
where  the  amount  of  additive  has  been  adjusted  to  balance  improvements  in  mechanical 
properties  and  UV  stability  against  loss  of  in  band  transmission. 

Given  the  nature  of  the  final  structure  which  is  expected  to  begin  with  a  fairly  rigid  base 
sheet  on  which  is  the  lattice  is  fabricated,  we  expect  that  base  will  be  either  HDPE  or 
Poly  IR2.  Filling  or  leveling  the  lattice  with  either  low  melting  PE  or  paraffin  have  been 
our  materials  of  choice. 

2.3  Fabrication  Procedure 

The  basic  substrate  fabrication  procedure  for  multiple  layers  involves  the  following  steps. 
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(1)  Fabrication  of  master  and  embossing  tool 

(2)  Microembossing 

(3)  Release  layer  deposition 

(4)  Metal  deposition  and  differential  stripping 

(5)  Interlayer  deposition  and  transfer  (and  align)  of  structures  from  donor  substrate 
(polycarbonate)  to  substrate  of  choice  (polyethylene) 

We  highlight  the  important  steps  in  the  discussions  below.  The  microembossing 
procedure  is  a  previously  developed  Polaroid/MicroContinuum  trade  secret. 

2.3.1  Master  Fabrication 

The  procedure  used  for  fabricating  the  embossing  tool  master  is  as  follows: 

(1)  Coat  photoresist  on  glass 

(2)  Expose  with  3-beam  interference  pattern  with  the  pitch  determined  by  beam  angles 

(3)  Modulate  exposure  and  develop  to  control  degree  of  hole  opening  and  depth 

(4)  Create  embossing  photoresist  master 

Figure  2.3. 1.1  A  and  2.3.1. IB  are  SEMs  of  a  master  fabricated  in  the  above  fashion. 

2.3.2  Microembossing 

Microembossing  is  a  technology  for  "stamping"  lithographically  defined  features  into  a 
substrate  of  choice.  The  advantage  is  that  it  offers  very  high  resolution  and  is  highly 
reproducible.  In  addition,  we  have  argued  that  it  is  one  of  the  few  choices  for  large  area, 
low  cost  fabrication  of  single  layer,  a  potentially  multiple  layer  3-D  photonic  crystal 
structures.  While  not  applicable  to  this  program,  it  should  be  pointed  out  that  the 
embossing/liftoff  procedures  we  demonstrate  for  micron  scale  features  has  also  been 
demonstrated  fori  00  nm  size  features. 

We  successfully  embossed  polyethylene  using  a  thermal  process,  and  CAP  and 
polycarbonate  using  a  chemical  process.  Both  of  these  processes  were  developed  using 
processes  that  are  proprietary  to  Polaroid  and  MicroContinuum. 

2.3.3  Release  Layer 

The  central  purpose  of  the  release  layer  is  to  allow  for  regions  of  the  metallic  overlayer  to 
be  precisely  and  repeatably  separated  from  the  underlying  substrate  leaving  behind  a  3-D 
periodic  structure.  We  performed  an  exhaustive  literature  search  and  many  trial  and  error 
attempts  to  determine  the  proper  release  agent.  The  best  results  were  obtained  using  a 
dilute,  aqueous  solution  of  cationic  detergent  that  coats  the  surface  with  a  very  thin  (sub 
micron)  residue  allowing  separation  of  the  metal.  Our  liftoff  experiments  were 
performed  with  Indium.  Indium  was  originally  chosen  because  it  has  a  low  vaporization 
temperature  and  seemed  to  be  most  suitable  for  differential  stripping  evaluation  studies. 
We  believe  the  use  of  detergent-based  release  agents  can  be  used  to  differentially  strip 
most  metals  using  the  techniques  described. 
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2.3.4  Differential  Stripping 


Differential  stripping  is  performed  by  pressing  the  donor  substrate  and  adhesive  layer 
between  2  rollers  and  then  pulling  them  apart.  We  tried  a  number  of  different  rolling 
pressures  and  adhesives  to  determine  the  proper  combinations.  We  found  that  regular 
scotch  tape  was  the  best  adhesive.  This  appeared  to  be  almost  pressure  insensitive. 

Figure  2.3 .4.1  illustrates  the  differential  stripping  process,  and  the  importance  of 
geometry  to  produce  the  required  effect.  In  particular,  the  process  works  because  the 
adhesive  does  not  come  into  contact  with  the  metal  in  the  cups,  and  hence  is  left  behind 
when  the  adhesive  is  removed.  Figure  2.3. 4.2 A  shows  a  surface  that  has  been 
differentially  stripped.  Figure  2.3.4.2B  shows  the  complimentary  liftoff  layer.  It  should 
be  pointed  out  that  this  process  was  extremely  reproducible  and  we  were  able  to 
demonstrate  excellent  differential  stripping  over  virtually  the  entire  sample  area  (~4  in2). 

The  importance  of  geometry  in  the  differential  stripping  process  cannot  be  over 
emphasized.  We  had  a  series  of  failures  in  the  stripping  process  that  were  eventually  tied 
to  a  faulty  master.  Figures  2.3. 4.3  A  and  2.3.4.3B  is  an  SEM  of  the  master  used  to 
fabricate  structures  that  did  not  adequately  strip.  The  embossing  master  shows  the 
features  that  will  be  embossed  in  a  substrate  using  the  corresponding  embossing  tool. 
Although  we  experienced  a  major  delay  associated  with  this  failure,  it  caused  us  to 
carefully  examine  the  stripping  process.  We  concluded  that  the  process  works  only  for 
"sharply"  defined  aspect  ratios.  Differential  stripping  works  for  the  hemispherical 
geometries  we  explored  because  less  material  is  deposited  along  the  ring  of  the 
hemisphere  allowing  it  to  fracture  in  the  region  of  interest,  and  thus  differentially  strip  as 
needed. 

2.3.5  Layer  Stacking 

Two  issues  dominate  the  procedures  that  can  be  considered  for  stacking  of  individual 
layers  to  form  a  3-D  structure.  First,  the  substrate  thickness  on  the  layers  to  be  stacked, 
such  as  those  shown  in  figure  2.3.4.2A  are  ~1  mil  (25  pm),  which  is  much  thicker  than  a 
lattice  constant.  And  second,  polyethylene  substrates  of  interest  do  not  have  good 
dimensional  stability.  This  severely  limits  the  ability  to  effectively  align  the  substrates. 

Considering  these  limitations,  we  developed  a  technique  for  stacking  structures  that  we 
refer  to  as  "direct  transfer".  The  procedure  is  as  follows: 

(1)  Begin  with  the  single  layer  structures  as  shown  in  figure  2.3.4.2A 

(2)  Fill  cups  with  acqueos  dispersion  of  polyethylene  polymer  coating  using  wire-wound 
coating  rod.  The  filling  is  controlled  by  the  solution  concentration  and  wire  guage  of 
the  coating  rod. 

(3)  Transfer  to  receiving  sheet. 

Figure  2. 3. 5.1  illustrates  the  procedure. 

Figure  2.3. 5.2  are  structures  that  has  been  transferred  to  a  receiving  sheet  using  the 
procedure  described.  Figure  2.3.5.3  is  a  cross  sectional  view  the  substrate  that  has  been 


planarized  as  described  in  step  (2)  above.  Figure  2.3.5.4A  and  2.3.5.4B  show  the  size 
over  which  structures  can  be  transferred.  Each  of  the  features  represents  a  "cup"  that  has 
been  transferred  from  the  donor  to  the  receiver  sheet. 

We  had  a  great  deal  of  success  using  this  procedure.  In  particular,  we  found  it  was 
relatively  easy  to  control  the  interlayer  thickness  as  well  as  to  reproducibly  transfer  large 
areas  from  the  donor  sheet  to  the  receiving  sheet. 

The  procedure  has  important  advantages  compared  to  other  approaches  we  considered. 
Most  importantly,  since  the  microembossed  structure  is  not  incorporated  into  final 
photonic  crystal  structure,  it  can  be  optimized  for  it's  microembossing,  mechanical,  and 
temperature  properties.  We  believe  that  the  transfer  procedure  is  the  most  important 
contribution  to  the  HIDE  program  by  the  SAIC  team. 

We  built  up  3-D  photonic  crystal  structures  using  this  procedure.  Figure  2.3. 5. 5  shows  a 
multilayer  stack  that  has  been  built  up  by  combining  the  steps  discussed  above.  To 
obtain  this  image,  we  cut  a  multilayer  structure  with  a  razor  blade  and  looked  for  regions 
that  had  been  cut  and  lifted  up.  Multiple  regions  can  be  seen  along  the  cut  line.  For 
example,  figure  2.3. 5. 6  shows  a  region  where  the  top  layer  has  separated  and  curled  back 
from  the  layers  below.  The  features  that  look  surprisingly  like  octopus  tentacles  are  the 
underside  of  the  upper  layer.  Figure  2. 3. 5. 7  is  a  high  magnification  image  of  a  single 
layer.  Notice  how  thin  the  interlayer  spacing  can  be  made. 

Figure  2.3. 5.8  shows  a  potential  difficulty  with  the  fabrication  procedure  described.  We 
were  concerned  that  the  stacking  process  would  crush  the  structures  below  the  top  layer. 
Figure  2.3. 5. 8  was  created  by  applying  the  same  pressure  and  heat  treatment  (both  of 
which  are  low)  used  in  building  up  the  structure  layer  by  layer.  The  hemispherical 
structures  are  flattened  somewhat.  We  explored  the  possibility  of  changing  the  filling 
material,  however  this  was  not  possible  in  view  of  the  time  and  budgetary  constraints. 

2.3.6  Alignment 

Alignment  of  the  layers  is  the  greatest  challenge  to  the  fabrication  procedure  described. 
The  initial  goal  was  to  create  zincblende  lattice  structures  by  stacking  individual  layers 
arranged  in  a  hexagonal  array.  This  approach  requires  6  layers  to  complete  a  unit  cell 
with  the  zincblende  structure  "grown"  in  the  [1 1 1]  direction.  Figure  2.3.6. 1  describes  and 
illustrates  the  stacking  procedure. 

We  explored  several  approaches  to  layer  alignment.  The  primary  difficulty  is  the  fact 
that  the  layer  and  substrate  must  be  exceptionally  rigid  to  allow  for  submicron  alignment 
over  wafer  size  substrates.  We  settled  on  a  procedure  that  involved  transfer  of  structures 
to  a  rigid  glass  substrate  and  aligning  subsequent  layers  to  software  generated  "fiducials". 
There  steps  were: 


(1)  Fuse  bottom  layer  to  glass  base. 


(2)  Align  next  layer  by  racking  focus  between  top  and  bottom  layers  and  using  software 
ficucials  and  features  (defects)  at  the  edge  of  the  pattern.  This  specific  step  is 
possible  because  of  the  many  reproducible  features  at  the  edge  of  the  pattern.  These 
features  are  the  same  for  all  layers  fabricated  from  the  same  master. 

(3)  Fuse  layers  and  repeat. 

(4)  Transfer  structure  off  of  rigid  glass  substrate  using  predeposited  release  layer. 

Figures  2.3.6.2A  and  2.3.6.2B  show  several  fiducials  that  are  marked  and  aligned  by  a 
software  image  processing  program. 

Maintaining  the  layer  dimensions  (i.e.,  not  stretching)  while  properly  holding  the  sample 
so  that  it  can  be  maneuvered  into  position  are  mutually  exclusive  and  we  developed 
several  generation  of  alignment  tools.  Figure  2.3.6.3  shows  the  final  setup  used.  This 
configuration  allowed  us  to  examine  the  diffraction  pattern  of  a  He-Ne  laser  impinging 
on  the  surface,  as  well  as  obtain  high  magnification  images  of  surface  features.  As 
mentioned  above,  these  features,  which  are  reproduced  on  each  layer,  were  used  as 
fiducials  for  layer  alignment.  Figure  2.3. 6.4  is  a  rendered  drawing  of  the  alignment  tool. 
The  tool  allowed  for  both  translational  and  rotational  alignment. 

We  performed  numerous  alignment  experiments  using  a  variety  of  techniques.  Our 
conclusion  is  that  while  rotational  alignment  is  possible,  translational  alignment  is 
extremely  difficult. 

3.0  Test  and  Measurement 

We  measured  both  single  layer  and  multilayer  structures.  The  single  layer  structures 
showed  general  agreement  between  predictions  and  measurements.  Figure  3.0.1  shows 
the  results  for  the  polyethylene  substrate  and  interlayer  (recall  the  interlayer  is  a 
polyethylene  aqueous  solution)  and  for  the  single  layer  hex  structures  (such  as  those 
shown  in  section  2.3.5).  The  high  angle  disagreement  between  experiment  and  theory  is 
likely  due  to  surface  roughness  and  the  fact  that  the  sample  is  partially  blocked  by  the 
sample  holder  at  high  angles. 

These  measurements  were  performed  using  a  CO2  laser  and  selecting  different 
wavelengths  within  the  available  9-1 1pm  band.  A  sketch  of  the  setup  is  shown  in  figure 
3.0.2. 

We  also  performed  measurements  of  the  multlayer  structures  (in  collaboration  with  SRI). 
These  results  were  inconclusive. 

4.0  Modeling  and  Simulation 

The  theoretical  effort  focused  on  formulating  a  semi-analytic  integral  equation  for  the 
computation  of  the  reflection  and  transmission  coefficients  of  a  periodic  single  or 
multilayer  array  of  dielectric  or  metallic  spheres  inside  a  host  medium,  which  may  be 
lossy.  The  results  of  this  effort  show  that  this  method  is  equivalent  to  that  of  applying  the 
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boundary  conditions  at  the  surface  of  the  spheres,  but  it  has  the  advantage  that  it  provides 
an  easy  transition  from  the  real  to  the  spectral  domain,  where  the  convergence  of  the 
computed  quantities  is  much  faster.  In  a  three  dimensional  lattice,  we  developed  two  for 
computing  the  dispersion  relation  of  the  photonic  crystal.  The  reflection  and 
transmission  coefficients  can  be  computed  for  a  single  or  multilayer  array  of  spheres 
embedded  in  a  slab.  Due  to  the  small  size  of  the  matrices  involved  in  the  computations 
and  the  fast  convergence,  there  is  a  small  requirement  for  computer  memory  and  time. 
The  attached  article,  which  has  not  been  submitted  for  publication  summarizes  the 
theoretical  results. 


Figure  2.1.1 


Figure  2.1.1:  Polyethylene  transmission  spectra 


Figures  2.3. 1.1 


Figures  2.3.1.1A  and  B:  Microembossing  masters 


Figure  2.3.4. 1 


Figures  2. 3.4.2 
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Figures  2. 3. 4. 3 


Figures  2.3.4.A  and  B:  Embossing  masters  used  to  emboss 
substrates  that  did  not  differentially  strip. 


Figures  2.3.5. 1:  Procedure  for  fabricating  multilayer 
structures 


a  donor  sheet  to  a  receiving  sheet  using  the  direct  transfer 
process 
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Figure  2.3.5.5:  Multilayer  structure 


Figure  2.3. 5. 6:  Lifted  off  top  layer  shows  how  thin 
the  interlayer  spacing  can  be  made 


Figure  23.5.1 :  High  magnification  of  the  top  layer 
illustrates  how  thin  the  interlayer  can  be  made 


Figure  2.3. 5. 8:  The  stacking  process  can  potentially 
crush  the  features  in  lower  layers 
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Figures  2.3. 6. 2 


Figures  2.3.6.2A  and  B:  Sample  fiducials  used  to  align 
layers 


Figures  2. 3. 6. 3 


Figures  2.3.6.  3:  Alignment  tool 
used  for  multilayer  alignment 


Figures  2. 3. 6.4 


Figures  23.6.4  :  Alignment  Tool  Detail 


Figures  3.0.1 
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Figures  3.0.2 


“  mini  BRDF”  (performed  at  NIST, 
Boulder) 


Reflection  and  Transmission  Coefficients  by  Periodic  Array 
of  Dielectric  Spheres  in  the  Spectral  Domain 


Abstract 

A  semi-analytic  integral  equation  approach  is  formulated  for  the  computation  of  the 
reflection  and  transmission  coefficients  of  a  periodic  single  or  multilayer  array  of  dielectric  or 
metallic  spheres  inside  a  host  medium,  which  may  be  lossy.  It  is  shown  that  this  method  is 
equivalent  to  that  of  applying  the  boundary  conditions  at  the  surface  of  the  spheres,  but  it  has  the 
advantage  that  it  provides  an  easy  transition  from  the  real  to  the  spectral  domain,  where  the 
convergence  of  the  computed  quantities  is  much  faster.  In  a  three  dimensional  lattice,  two 
methods  are  presented  for  computing  the  dispersion  relation  of  the  photonic  crystal.  The 
reflection  and  transmission  coefficients  are  computed  for  a  single  or  multilayer  array  of  spheres 
embedded  in  a  slab.  Due  to  the  small  size  of  the  matrices  involved  in  the  computations  and  the 
fast  convergence,  there  is  a  small  requirement  for  computer  memory  and  time. 

I.  Introduction 

In  recent  years  there  has  been  significant  progress  in  the  theoretical  and  experimental 
investigation  of  photonic  band  gap  crystals.  In  most  of  the  work  that  has  been  done,  simple 
shapes  for  the  dielectric  or  metallic  materials  were  chosen  such  as  cylinders  with  circular  or 
rectangular  cross-section  for  the  two-dimensional  crystals  and  spheres  or  rods  with  rectangular 
cross-section  for  the  three-dimensional  crystals.  The  scattering  of  light  by  a  single  layer  periodic 
array  of  spheres  was  treated  by  K.  Ohtaka  et.  al.  [1-4],  where  the  analogy  between  low-energy 
electron  scattering  and  the  light  scattering  by  the  dielectric  spheres  was  demonstrated  and  the 
connection  was  made  between  the  reflection  of  light  and  its  integrated  density  of  states.  A 
formalism  for  the  same  problem  was  developed  by  A.  Modinos  [5]  by  the  direct  application  of 
the  boundary  conditions  at  the  surface  of  the  spheres.  This  formalism  was  applied  in  the 
computation  of  the  reflectivity  of  a  single  layer  array  by  N.  Stefanou  et.  al.  [6]  and  in  the 
computation  of  the  transmittance  of  a  multilayer  array  of  the  photonic  band  structure  of  a  three- 
dimensional  array  by  A.  Modinos  et.  al.  [7]. 

Various  techniques  have  been  developed  for  the  computation  of  the  dispersion  relation  of  the 
photonic  crystal  such  as  the  plane- wave  expansion  method  [8-12]  and  the  transfer  matrix  method 
[13-16].  The  former  method  has  been  applied  for  scalar  and  vector  fields.  In  the  vector  case,  the 
fields  are  discontinuous  at  the  boundaries  between  different  dielectrics.  Consequently,  the 
truncation  of  the  Fourier  series  in  the  spectral  domain  contains  high  spectral  components,  so  that 
convergence  becomes  a  serious  problem  for  the  vector  fields.  For  example,  R.D.  Meade  et.  al. 
[11]  used  106  plane  waves  to  obtain  reliable  numerical  results.  Convergence  is  improved  if  the 
sharp  edges  are  smoothed  to  suppress  the  high  spatial  frequencies  in  the  Fourier  series.  This  was 
done  by  H.S.  Suziier  et.  al.  [12]  to  remove  the  discontinuity  of  the  fields.  Even  so,  the  plane 
wave  method  is  not  sufficient  for  more  elaborate  photonic  crystals  such  as  crystals  with  defects 
or  containing  metal.  The  transfer  matrix  method  is  applied  in  real  space  and  is  fairly  easy  to 
implement.  It  provides  for  the  computation  of  either  the  dispersion  relation  or  the  reflection  and 
transmission  coefficients  for  crystals  of  finite  size.  Due  to  the  large  number  of  mesh  points 
needed  for  accurate  computations,  the  transfer  matrix  method  may  lead  to  numerically  unstable 
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solutions.  Another  method  that  was  developed  to  avoid  this  problem  is  the  R-matrix  propagator 
method  [17].  It  was  applied  to  two-dimensional  structures  to  compute  both  the  band  structure 
and  the  transmission  in  photonic  media. 

In  this  paper,  the  scattered  field  and,  therefore,  the  reflection  and  transmission  coefficients 
for  a  single  layer  periodic  array  of  dielectric  spheres  embedded  in  a  host  medium  are  computed 
with  the  help  of  an  integral  equation  which  is  known  in  electromagnetic  theory  as  the  Ewald- 
Oseen  extinction  theorem  [18].  The  unknown  quantities  which  are  determined  by  using  this 
theorem  are  the  electric  and  magnetic  currents  at  the  surface  of  the  spheres.  Then  the  fields 
inside  and  outside  the  spheres  can  be  computed  from  these  surface  currents.  In  the  integral 
equation  approach  presented  in  this  paper,  either  the  electric  field  by  itself  or  the  magnetic  field 
by  itself  is  sufficient  for  the  whole  computation  of  the  scattered  field.  This  is  not  the  case  when 
the  boundary  conditions  are  applied  where  both  the  electric  and  magnetic  tangential  components 
to  the  surface  of  the  spheres  are  needed  for  the  computation.  The  integral  equation  method  is 
shown  to  be  equivalent  to  that  by  A.  Modinos  [5],  but  it  leads  naturally  into  formulating  the 
problem  from  the  real  into  the  spectral  domain.  In  a  lossless  host  medium,  the  main  advantage 
of  the  spectral  domain  computation  is  that  it  provides  faster  convergence.  Thus,  in  a  two- 

dimensional  array,  the  terms  to  be  summed  in  the  function  Z]  ™  (cf.  Section  IV)  are  inversely 
proportional  to  the  distance  of  the  lattice  in  the  real  domain,  while  they  are  inversely  proportional 
to  the  cubic  power  of  the  distance  in  the  reciprocal  lattice  (spectral  domain).  In  a  three 
dimensional  array  it  is  even  better,  since  the  terms  to  be  summed  are  still  inversely  proportional 
to  the  distance  of  the  lattice  in  the  real  domain,  but  they  are  inversely  proportional  to  the  fourth 
power  of  the  distance  in  the  reciprocal  lattice.  In  a  lossy  host  medium,  if  the  loss  is  fairly  large, 
the  computation  becomes  faster  in  the  real  domain  due  to  the  exponentially  decaying  terms 
present  in  the  expressions  for  the  fields. 

In  the  integral  equation  approach  and  in  the  approach  by  A.  Modinos  [5],  the  fields  are 
expanded  in  vector  spherical  harmonics  with  (l,m)  eigennumbers.  In  Table  I  in  the  paper  by  N. 
Stefanou  [6]  it  is  shown  that  a  maximum  value  of  lmax=5  is  sufficient  for  accurate  results.  The 
dimension  of  the  matrices  involved  in  the  computations  is  21max(lmax+2)*2Imax(lmax+2).  Thus, 
with  lmax=5,  the  dimension  of  the  matrices  is  70x70  which  is  modest  compared  to  the  plane  wave 
method.  Consequently,  the  scattering  and  transfer  matrices,  at  frequencies  where  only  a  few 
propagating  modes  are  present,  are  small  in  size  and  they  should  not  lead  to  unstable  solutions 
when  multilayer  arrays  are  considered.  Another  advantage  of  the  integral  equation  approach, 
since  it  is  equivalent  to  that  of  applying  the  boundary  conditions,  is  that  the  discontinuity  of  the 
fields  at  the  sharp  boundaries  (surfaces  of  the  spheres)  is  built  into  the  solution,  which  is  not  the 
case  in  the  plane-wave  method. 

In  a  three-dimensional  lattice,  two  approaches  are  provided  for  the  computation  of  the 
dispersion  relation  of  the  photonic  crystal  (cf.  section  V),  either  by  setting  the  determinant  of  a 
linear  homogeneous  algebraic  system  equal  to  zero  or  by  applying  Bloch’s  theorem  as  it  was 
done  in  the  paper  by  N.  Stefanou  et.  al.  [6].  Since  the  size  of  the  matrices  involved  is  small,  the 
requirement  in  computer  memory  is  small.  Also,  from  the  discussion  above,  it  follows  that  the 
computation  in  the  spectral  domain  should  be  fast  and,  therefore,  the  requirement  in  computer 
time  is  small.  The  main  disadvantage  of  the  method  presented  in  this  paper  is  that  it  is  confined 
to  non-overlapping  dielectric  or  metallic  spheres.  But  the  method  can  be  extended  to  other 
shapes  of  scatterers  as  long  as  a  complete  set  of  functions  exists  to  expand  the  surface  electric 
and  magnetic  currents.  Thus,  the  method  can  be  extended  to  cylinders  of  circular  or  rectangular 


cross  section  and  to  multilayer  periodic  structures  where  each  layer  consists  of  parallel  rods  of 
circular  or  rectangular  cross  section. 

II.  Integral  Equation  Approach  for  a  Single  Sphere 

A  dielectric  sphere  is  centered  at  the  origin  of  the  coordinate  system  and  has  conductivity  g3 
and  dielectric  constant  s3  (in  Mks  units).  It  will  be  referred  as  medium  3.  Medium  2  outside  the 
sphere  has  conductivity  g2  and  dielectric  constant  82.  We  define  the  complex  wavenumbers  Oj, 
i=2,3  as  follows: 

<*?  =  j®Po(gi  +j®eiX  (1) 

where  po  is  the  magnetic  permeability  of  vacuum,  co  is  the  frequency  and  j  is  the  imaginary  unit. 
Everywhere,  in  the  following,  the  time-dependence  is  assumed  to  be  e1®1.  In  terms  of  Cj,  the 
outward  radial  Green’s  functions  in  the  two  media  are  defined  by  the  relations 


Gj(r) 


47ir 


-e-°ir. 


where  r=|r|  and  i=2,3. 

At  the  surface  of  the  sphere. 


we  define  the  electric  and  magnetic  surface  currents 
J<e>(S)  =  [erxH3(r)]r=R, 

J(m)(S)  =  -[er  x  E3(r)]r=R, 


(2) 


(3a) 

(3b) 


Here,  E3(r),  H3(r)  are  the  electric  and  magnetic  fields  inside  the  sphere,  S  is  a  point  on  the 
surface  SR  of  the  sphere,  determined  by  the  spherical  coordinates  (R,  0,  9),  R  is  the  radius  of  the 
sphere  and  er  is  the  outward  unit  vector  normal  to  the  surface  of  the  sphere. 

In  the  following  we  shall  use  the  vector  spherical  harmonics  [19]  as  the  complete  set  to 
express  the  angular  dependence  of  vectors  on  Qr  =(0,  9)-  Their  definition  and  the  relations  they 

satisfy  are  given  in  Appendix  I.  Since  YIR(Qr)is  normal  to  er,  as  well  as  erx  Y,™(Qr),  the  most 
general  expression  of  the  electric  and  magnetic  surface  currents  is  as  follows: 


j'"(s) = dCviSca)  -  Ce,  x  v,:]. 

Im 


(4a) 


J<M)(S)  =  £lvi!nEX(nt)- jV<mMle,  X  Y,™(Qr)]>  <4b> 

lm 

where  the  summation  extends  over  1=1, 2, 3,...,  m=-l, -1+1,..., 1-1,1.  The  four  sets  of  unknown 

coefficients  l{^\  l{m\  V,^,  be  determined  by  solving  an  integral  equation  and 

imposing  one  of  the  conditions  in  Eq.  (3a)  or  Eq.  (3b). 

Let  us  define  the  electric  and  magnetic  vector  potentials  by  means  of  the  surface  integrals  of 
the  electric  and  magnetic  surface  currents,  i.e.. 


(5) 


A[h)(r)=  jGi(r-r')J(h)(S')dS', 

Sr 

where  h=e,m  and  i=2,3.  In  accordance  with  the  paper  by  E.  Wolf  [18],  we  define  the  electric  and 
magnetic  fields  for  i=2,3  namely, 


E-J)(r)  =  -^V  x  V  x  A[e)(r)  +  V  x  A[m)(r), 

CTi 

HjJ)(r)  =  -  V  x  V  x  Alm)  (r>  -  V  x  Aie>(r)’ 


Since 


V2A-h)(r)-C7?A[h)(r)  =  0, 

for  all  r  except  for  r  on  Sr,  an  equivalent  expression  for  these  fields  is: 


E[J)(r)  =  jwpo 


H[J)(r)  =  -p- 


Aje)  (r)  — y  V(V  •  A-e>  (r) 


+  Vx  A-m)(r), 


Jo^O 


M)  (r)  — -fV(v  •  A[m)(r)) 


-  Vx  Aje)(r), 


(6a) 

(6b) 


(7) 


(8a) 


(8b) 


The  superscript  (J)  indicates  that  EjJ)(r),H[J)(r)  depend  on  J(e)(S),  J(m)(S)  and  therefore,  they 

depend  on  the  unknown  sets  of  coefficients  l{^  \  l}„ ,  V,^,  V® . 

First,  we  confine  ourselves  in  the  region  inside  the  sphere.  Substitution  of  the  Green’s 
function  as  given  by  Eq.  (II 1)  in  Appendix  I  and  the  surface  currents,  as  given  by  Eqs.  (4a),  4b), 

and  a  rather  lengthy  but  straightforward  computation  of  E-^(r),  H^(r)  in  Eqs.  (6a)  and  (6b) 
leads  to  the  following  relations: 


lm 

(9a) 

HjJ,(r)  -  jJ-Z  fe-'ESlrV)  +  fffl"”ESn|(r)], 

lm 

(9b) 

E!S'"'(r)  =  f,(cir)Y,;(Qr)0(R  -  r). 

(10a) 

ES”'(r)  =  ]i-Vx(fl(0lr)Yi;(Qr))©(R  -  r), 

(10b) 

where  i=2,3  and 


The  function  @(x)  is  defined  as  follows: 


Also: 


@(x)  = 


1  x  >  0, 
0  x  <  0 . 


r-(M,in)  _  2 

ri,lm  _  zi 


“HO 


(11) 


(12a) 


/r(E,in)  _  2 
ri,lm  “  zi 


where  Zj=a-,R  and 


l[zg,(z)J  =g;  (z)  +  {gi(z)  =  ~2i+rgi+i(z)  “  2iTTgi-i(z)- 


(12b) 


(13) 


The  prime  indicates  differentiation  with  respect  to  z.  Use  was  made  of  Eq.  (13)  in  deriving  Eqs. 
(12a)  and  (12b). 

The  fields  inside  the  sphere  are  given  by  E3(r),  H3(r).  Notice  that  these  fields  are  not 
independent  of  each  other,  since  they  satisfy  Maxwell’s  equations.  If  one  of  the  two  fields  is 
known,  the  other  can  be  computed  from  the  curl  of  the  known  field.  Next,  one  of  the  conditions, 
given  by  Eqs.  (3a)  and  (3b)  will  be  imposed  on  the  field  at  the  surface  of  the  sphere.  Applying 
the  identities  (I5a)-(I5c)  in  Appendix  I,  on  Eqs.  (9a)  and  (9b)  we  obtain  the  relations  (we  omit 
the  superscript  (J)  for  the  fields  inside  the  sphere). 


lm  *" 

[e,  x  H,(r)U  =  f&n)irEf,(^)l  YS(£2r)  -  *  Y^fi,) 


lm 


(14a) 

,  (14b) 


where  use  was  made  of  the  identity. 

7 [zf, (z)|  s  f,'(z)  +  4f((z)  =  2i+Tfi+i(z)+ 2i+r fi-i (z>’ 
Condition  (3  a)  provides  the  relations 

Cin,Thf,(Z3)U^C\ 


°3 


(15) 


(16a) 


5 


(16b) 


E(E,in)f  f  ,_®^i(E) 
fc3,lm  ~  CT,  l\m  ■■ 


while  condition  (3b)  gives  the  relations: 


r(M,in)f  ,  \  y(M) 

t3,lm  rlW  vlm  ’ 


(17a) 

(17b) 


where  E^"1*,  are  given  by  Eqs.  (12a)  and  (12b)  for  i=3.  If  use  is  made  of  the  identity 


^-  +  ^[zgl  (z)J  f|  (z)  =  i[zfi  (z)J  g ,  (z), 

either  set  of  Eqs.  (16a),  (16b)  or  Eqs.  (17a),  (17b)  gives  the  same  relations,  namely 


y(M) 

vlm 

fl  (z3  ) 

c>H0  i(M) 
ct3  lm 

^-[z3f,(z3)] 

^-[z3f,(z3)] 

[(E) 

f|(z3)  ' 

(19a) 


(19b) 


Thus,  the  four  sets  of  unknown  coefficients  have  been  reduced  to  two  sets.  We  introduce  the  two 
independent  sets  A[„  \  A{^  by  means  of  the  relations 


C=^3f3fe)fALM> 


(20a) 

(20b) 

(20c) 

(20d) 


Substituting  these  expressions  into  Eqs.  (12a),  (12b).  for  i=3,  and  making  use  of  Eq.  (18),  we 
obtain  the  relations 


-(M,in)  _  a(M) 


(21a) 


6 


w  B 

II 

'Hz 

5-1 

XT 

(21b) 

• 

On  the  other  hand,  for  i=2,  we  obtain  the  relations 

/r(M,in)  _  •  (M)a(M) 
t2,lm  ~  al  Alm  ’ 

(22a) 

• 

r-(E,in)  _  (E)a(E) 

^2,lm  ~  al  Alm  ’ 

(22b) 

where 


„(M)  __7 

a,  -  z2 


fl  (Z3)[Z2§I  (Z2>1  "  §l(Z2  )IZ3fl  (Z3)l 


(23a) 


aiE)  =  -z, 


fl(z3)[z2gl(Z2)J-(^-ygl(z2)[Z3fl(Z3l 


(23b) 


Therefore,  in  terms  of  the  unknown  coefficients  the  fields  inside  the  sphere  are 

given  by  the  usual  expressions 


E,(r)  =  iElAir'Effi-’fr)  +  A^'W]. 

Im 


(24a) 


H3(r) 


CT3 


(E)F(M,in) 
jtono  L"  3,1m 
Im 


A^E™(r)  +  ArEf 


(24b) 


The  two  sets  of  unknown  coefficients  a{^  ,  A  are  determined  by  an  integral  equation, 

which  is  the  Ewald-Oseen  Extinction  Theorem  for  the  electric  field  (or  for  the  magnetic  field),  as 
stated  in  the  paper  by  E.  Wolf  [18].  The  theorem  states  that  for  all  r  inside  the  surface  S  (in  our 
case  inside  the  sphere)  the  following  integral  equation  is  satisfied: 


E^J)(r)  +  Ef  (r)  =  0,  (25) 

where  E^OO  is  given  by  Eq.  (6a)  or  Eq.  (8a)  in  general  and  by  Eqs.  (9a),  (22a),  (22b),  in 
particular,  for  the  case  of  the  sphere.  Therefore,  for  the  sphere,  E^1  (r)  is  equal  to 

E<2’>(r)  =  -jX[a]M'A!“»E«“r>(r)  +  a<E'A^»ES;'(r)j.  (26) 

Im 

E(,S)(r)  is  the  source  (or  driver)  electric  field  that  originates  outside  the  surface  S,  i.e.,  in 
medium  2. 
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Similarly,  the  unknown  coefficients  can  be  determined  by  the  Ewald-Oseen 

Extinction  Theorem  for  the  magnetic  field  [18],  which  states  that  for  all  r  inside  the  surface  S, 
i.e.,  the  sphere,  the  following  integral  equation  is  satisfied: 

H^J)(r)  +  H(2S)(r)  =  0,  (27) 

where  H^(r)  is  given  by  Eq.  (6b)  or  Eq.  (8b)  in  general  and  by  Eqs.  (9b),  (22a),  (22b),  in 
particular,  for  the  sphere.  Then,  for  the  sphere,  H^(r)  is  equal  to 

H<»(r)  =  -^S[a!EIA!?Effi"»(r)  +  a|M>  A<“'E'25*>(r)].  (28) 

l,m 


H^(r)  is  the  source  magnetic  field  that  originates  outside  the  surface  S. 

If  the  source  field  is  a  free  field  in  medium  2,  then,  in  the  plane  wave  representation,  it  is  given 
by  the  expression 

E<s)(r)  =  e~jk p  p |e(ute+) (0)e~Y 2 z  +  E(kTE_)(0)eY2Z|ez  x  kp) 

+  (E(k™+)(0)e“Y2Z  +  E(k™-)(0)eY2Z)kp  (29a) 

-  (Ek™+)  (0)e_Y  2  z  -  Ek™_)  (0)eY  2  z  ez " , 


where 


H(2S) 


(r)  =  ^2_e 
v  ;  jwno 


-jkpP 


-  (E(kTE+) (0)e“Y  2 z  -  EkTE_)  (0)eY 2 z )^-! 


+  (E(k™+)(0)e‘Y2Z  -E<™-)(0)eY2Z)^-(ez  xkp) 
+  (Er)(0)e-Y2Z  +  E<tTE->(0)eY2Z)^ez', 


(29b) 


Y2  = 


(30) 


Here,  EkTE±)(0),  E(kTM±)(0)  are  the  amplitudes  of  the  electric  field  at  z=0,  moving  in  the  forward 
(+)  and  backward  (-)  directions,  for  the  transverse  electric  (TE)  and  transverse  magnetic  (TM) 
plane  waves,  respectively.  Also,  kp  is  the  projection  of  the  wavevector  on  the  xy-plane,  kp  is  its 

magnitude  and  kp  is  the  unit  vector  in  the  direction  of  kp.  If  (0k,  (pk)  provide  the  direction  of  the 
plane  wavefront  with  constant  phase,  in  spherical  coordinates,  then 


kp  =  ex  cos  (pk  +  ey  sin  cpk , 


(31a) 


kp  =  ksin9k, 


(31b) 


y2i=kcosek,  (31c) 

where  y2i,  is  the  imaginary  part  of  y2  and  k  is  to  be  determined.  If  o2  is  expressed  by  the 
complex  index  of  refraction,  i.e., 

ct2  =(k2+jn2)t,  (32) 


where  c  is  the  velocity  of  light  in  vacuum,  then,  k  can  be  computed  with  the  help  of  Eqs.  (39), 
(31b),  (31c)  and  kp  becomes  equal  to 


=  Vn2_k 


2f  sin9k 


1  +  . 


1  + 


2k 2^2 


^  n  2  )cos®k 


Yi 


(33) 


For  simplicity,  we  assumed  that  k2<n2  which  is  usually  the  case. 

If  use  is  made  of  Eq.  (II  19)  in  Appendix  II  and  of  the  orthonormality  of  the  scalar  spherical 
harmonics,  then,  the  source  field,  in  the  spherical  wave  representation  becomes  equal  to 

Ef(r)  =  jE[Ar,fi92r)VlS(a)+  a£S) +Vx (f,(o2r)Y{;(Clt))|  ,  (34a) 

lm 

h2SV)  -  j^ZlAES>fi(CT2rW(n, )+ A«  X  (f,(02r)v;(nt))J ,  (34b) 

lm 

where 


and 


*{c!™kE*»(0)  +  (-l)'-“'E're->(0)] 

+ (f  Jfcj™1 1  [e(™*>(0) + (_iy— ■  E«™ ->(0)]}. 

a  (ES)  _  a2  (  i\l-m+l 

Alm  “  r  ,1/  Y2  V  L) 

[|(I+1)1  2  72 

*  {c}™’  [E(kTE+)(9)  +  (-l)1_m+lE(kTE_)(0)] 

-  C^E)[E(k™+)(0)  +  (_iy-^‘E<™  _)(0)]} , 


(35a) 


(35b) 
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p(TE) 

^lm 


=  j 


9k 


(36b) 


ar=i[(l-mXl  +  m  +  l)l^5  (37a) 

pr=i[(l  +  mXl-m  +  l)J^.  (37b) 

The  complex  scalar  spherical  harmonics  Ylm(^-,cpk)  have  been  defined  in  Appendix  II,  Eq.  (II 

19).  When  medium  2  is  non-absorptive,  i.e.,  02  is  purely  imaginary,  they  reduce  to  the  usual 
scalar  spherical  harmonics. 

Substitution  of  the  source  electric  field  in  the  spherical  wave  representation,  into  Eq.  (25) 
determines  the  unknown  coefficients  A® ,  namely, 


a(M)  I  a  (MS) 
'Nm  (M)  1m  ’ 
ai 


(38a) 


a(E)  -  I  A.(ES) 
'Mm  (E)  ’ 

“I 


(38b) 


where  use  was  made  of  the  orthonormality  property  of  the  vector  spherical  harmonics.  These  are 
the  same  expressions  obtained  by  the  application  of  the  boundary  conditions  at  r=R,  i.e., 
continuity  of  the  tangential  components  of  the  electric  and  magnetic  field.  Notice  that  the 
application  of  the  Ewald-Oseen  extinction  theorem  led  to  the  determination  of  the  unknown 

coefficients  a[^,  a|^  without  any  reference  to  the  fields  outside  the  sphere.  Also,  once  the 

electric  and  magnetic  surface  currents  were  expanded  in  terms  of  a  complete  set,  one  of  the 
boundary  conditions,  i.e.,  either  Eq.  (3a)  or  Eq.  (3b)  was  indispensable  in  reducing  the  four  sets 

of  unknown  coefficients  V^,  V,^\  ijj^,  I®  into  two  sets.  The  same  result,  i.e.,  Eqs.  (38a), 
(38b)  would  have  been  obtained  if  the  Ewald-Oseen  extinction  theorem  for  the  magnetic  field 
was  applied. 

Therefore,  in  general,  it  is  possible  to  determine  the  fields  in  a  dielectric  medium  3  which  lies 
inside  medium  2  with  their  boundary  being  the  surface  S  and  the  source  field  originating  outside 
medium  3.  For  this  purpose,  the  electric  and  magnetic  surface  currents  on  S  are  expanded  in 
terms  of  a  complete  set  of  functions  with  four  unknown  sets  of  coefficients.  If  use  is  made  only 
of  the  electric  field  inside  S,  then  the  boundary  condition,  Eq.  (36),  is  imposed  which  provides 
two  linear  algebraic  sets  of  equations  for  the  four  unknown  sets  of  coefficients,  by  making  use  of 
the  orthonormality  property  of  the  complete  set.  Two  more  linear  algebraic  sets  of  equations  are 
provided  by  the  integral  equation  (25).  These  four  sets  of  linear  algebraic  equations  with  four 
unknown  sets  of  coefficients  determine  the  electric  field  inside  S.  The  magnetic  field  inside  S  is 
computed  from  the  curl  of  the  electric  field.  On  the  other  hand,  if  use  is  made  only  of  the 
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magnetic  field  inside  S,  then  the  boundary  condition,  Eq.  (3a),  together  with  the  integral  equation 
(27)  determine  the  magnetic  field  inside  S.  The  electric  field  is  computed  from  the  curl  of  the 
magnetic  field. 

Next,  we  consider  the  scattered  field  outside  the  sphere  S.  Substitution  of  the  Green’s 
function,  i.e.,  Eq.  (Ill)  and  the  surface  currents,  i.e.,  Eqs.  (4a),  (4b),  into  Eqs.  (6a),  (6b)  for  i=2, 
leads  to  the  relations 


Ef  >(r)  =  jEkrr’E^r’W  +  ^hT’ESTV)],  (39a) 

1m 

Hf'Or)  -  +  ^"’E^rV)].  (39b) 

lm 

where 

E(ioutV)  =  g,(a2r)Y^(Qr)0(r  -  R),  (40a) 


ES7V)  -  j^V  x  (g,(cr2r)Y“(nt))s>(r  -  R), 

(40b) 

fa0""  =-4[“ f|(z2)C 

(41a) 

(41b) 

We  have  renamed 

the  fields  to 

indicate  that  they  are  the  scattered  fields. 

Replacing 

ij^),  I®,  yff  from  Eqs.  (20a)  -  (20d)  into  the  relations  above,  we  obtain 

r-(M,out)  _o(M)a(M) 
t2,\m  ~Pl  Alm  ’ 

(42a) 

^(E,out)  _  o(E)  .(E) 

C2,[m  ~  Pl  Alm  ’ 

(42b) 

where 

Pr=z2 

fl  (z3 )tz2fl  (z2  )J  -  fl  (z2  )tz3fl  (z3  )1 

9 

(43a) 

P{E)=Z3  fl(z3)[z2f,(z2)|  -(ij)2 f1(z2)[z3f1(z3)|  . 

(43b) 

The  total  field  outside  the  sphere  is  equal  to 


E2(r)  =  E<2s)(r)  +  Ef)(r). 

(44a) 

H2(r)  =  H<s)(r)  +  Hf)(r). 

(44b) 

When  the  source  field  is  a  free  field,  as  given  by  Eqs.  (29a),  (29b),  then,  with  the  help  of  Eqs. 
(38a),  (38b),  (42a),  (42b),  the  scattered  field  due  to  the  surface  currents  becomes  equal  to 


Ef  >  (r)  =  |BS><»r> «  +  b£)E^U°  (r)| 

i,m 

tt(sc),x_  ct2  vfR^E^F(M’outVr'l  + R^F^E’0UtVrJ 

H2  (r)_  2- Pirn  L2,lm  W  +  Blm  L2,lm  WJ 

u  l,m 


(45a) 

(45b) 


where 


n(M) 

R(M)  _  P| _  *  (MS) 

0lm  (M)  lrn  ’ 


(46a) 


r(E)  _  pi 

Dlm  (E)  ■r^lm 


(E) 


a<!S). 


(46b) 


These  are  the  same  expressions  obtained  by  applying  the  boundary  conditions  at  the  surface  of 
the  sphere.  Therefore,  the  integral  equation  approach  provides  the  correct  solution  inside  and 
outside  the  sphere.  In  the  general  case  of  any  surface  S,  the  difficulty  of  the  formalism  presented 
lies  in  the  choice  of  the  complete  set  to  expand  the  electric  and  magnetic  surface  currents.  In  the 
case  of  the  sphere,  the  obvious  choice  is  the  vector  spherical  harmonics,  which  led  to  a  trivial 
solution  for  the  linear  algebraic  set  of  equations.  As  we  will  see  in  the  next  section,  this  is  not 
the  case  for  a  two-dimensional  periodic  array  of  dielectric  spheres  embedded  in  a  host  medium. 

III.  Integral  Equation  Approach  for  an  Array  of  Spheres 

A  single  layer  periodic  array  of  dielectric  spheres  with  conductivity  gs  and  dielectric  constant 
S3  (medium  3)  is  embedded  in  the  host  medium  2.  The  unit  call  in  the  lattice  is  determined  by 
the  translation  vectors  ai,  a2  so  that  the  points  in  the  lattice  are  given  by  the  relation 


apq=pa,+qa2,  (47) 

where  p,  q  are  integers.  The  fundamental  vectors  of  the  reciprocal  lattice  are: 


B,  =27t  "2XC*  v, 

1  »l(a2xcz) 


B2  =  271  *z*a|  r  , 

1  al  (a2  xez ) 


and  the  points  in  the  reciprocal  lattice  are  given  by  the  relation 


(48a) 

(48b) 
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(49) 


bpq  =  PB1  +9B2- 

The  dielectric  spheres  are  centered  at  apq  and  they  do  not  overlap.  At  the  surface  Ssph  of  each 
sphere  there  is  an  electric  and  magnetic  surface  current.  Due  to  the  periodicity  of  the  system 
Floquet’s  theorem  applies,  and  therefore,  the  electric  and  magnetic  surface  current  due  to  all  the 
spheres,  centered  at  apq,  are  equal  to 


J"V)  =  Ze-*' *" ufr  -  ap„K6  («,„  )  -  jC%,  -  *1"  («,„  )1  (50a) 

pq  lm 

j(“>(I.)  =  2e'JP“Nu|r-aM|ESspJ*f;[vS>Y|S(nrM)-jV|<inMle,pq»<Y^(nrpq)]  (50b) 

pq  lm 


Here,  p  is  the  phase  vector  mentioned  above,  Qrpq  =  (0pq,9pq)>  where  (rpq,  0pq,  <ppq)  are 
the  spherical  coordinates  of  the  vector  rpq  =  r  -  apq ,  and 


pq 


1 

0 


if !  r  "apq  HR- 
otherwise 


(51) 


The  surface  current  on  each  sphere  satisfies  the  conditions  given  by  Eqs.  (3a)  or  (3b).  Using 
one  of  these  conditions,  and  following  the  same  procedure  as  in  Section  II,  we  conclude  that  Eqs. 
(20a)  -  (20d)  are  valid,  as  well  as,  Eqs.  (21a),  (21b)  and  Eqs.  (22a),  (22b).  Therefore,  the  fields 
inside  the  sphere  centered  at  apq,  i.e.,  when  |r-apq|<R,  is  given  by  Eqs.  (24a),  (24b),  where  r 
should  be  replaced  by  rpq.  On  the  other  hand,  when  |r-apq|<R  in  the  computation  of 

E^J)(r),  H(2J)(r),  the  field  due  to  the  sphere  centered  at  apq  is  computed  inside  that  sphere,  while 
for  all  the  other  spheres,  the  field  due  to  them  is  computed  outside  these  spheres.  Taking  into 
account  Eqs.  (42a),  (42b),  we  conclude  that 


(J)rrt  -  _;o~jP'aPci 


E^(r)  =  -je 

f 


lm 


(M) 

‘■lm 


a!M,Effi",4)-PiM,Z'e"Pip',E^r»(rp,-ap,,) 


\ 


V 


pq 


(52a) 


J 


+  A2>  ajE’E<“>(rpq)-(JiE>Z  e 


(E)^  4-jP'Vq'F(E.out) 

p'q 


E(2Elm°ut)(rpq-apV) 
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HV  >  (r)  =  ^2_e  ”J|J  ap<1 


X  A® 


+  ALM’  a!M)E®:'(rp,)-p!M»X,e-iN'’,‘''E®r(rpq  -5pV) 


(52b) 


where  ap-q.  =  ap.q.  -  apq .  The  prime  in  the  summations  above  means  that  they  extend  over  all 
spheres  except  for  the  one  centered  at  apq.  The  (in)  and  (out)  fields  in  the  relations  above  have 
already  been  defined  in  the  previous  section. 

Multiplication  of  Eqs.  (52a),  (52b)  by  (rpq  )  and  integration  over  the  volume  Vpq  of 

the  sphere  centered  at  apq  yields  the  following  relations: 


*ECi  ^nr8„.5mm.-a{i,.m.(p)  B|«>-ngrm.(P)B®  , 

I'm'  Pi  J  I 


(53a) 


J  HjJ>(r)  •  E^"**  (rpq)d3x[ 


G2  apq 


Tc,  -nffi, w(P)B!“’+  Sir8ll.8mm.-o®l.ln.(P)  b] 

i'm'  L  Pi 


(53b) 


where 


<i'm'(P)  =  rEe-jN^UIW;lm(ap,q,) 
Si  p'q' 


(54c) 


n£m,'  (P) =  r  E <riPV'1' u  Gw  )• 

SI  p'q' 


(54d) 
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(55a) 


U»W  fE<2Mr’(r-r0).E«S>  (r)d3x, 
Vo 


V,mJ.m.(r0)  =  jESS“»(r-r0)-B®i»  (r)d3x 
v0 


(55b) 


Ci  =  flfi(^2r)l2  f2dr 


[Z2fl-I (z2 )f| (Z2 )  -  z*2fl-l (2*2 )f| (Z2) ] • 

2 


(56) 


Here,  we  have  replaced  the  unknown  coefficients  by  =  Pim  )Aim  ^ 

=  PjE^Aj^>  since  the  latter  provide  the  scattering  properties  of  the  array  in  which  we  are 
interested.  Also,  8mn  is  the  Kronecker  delta,  i.e.,  Smn=l,  when  m=n  and  zero  otherwise. 

In  the  presence  of  the  lattice,  the  source  field’s  wavevector  is  not  limited,  in  general,  to  the 
wavevector  kp  of  the  incident  field.  Free  fields  with  wavevectors  Kpq=kp+bpq  where  bpq  is  a 
point  in  the  reciprocal  lattice  (Eq.  (49)),  could  also  drive  the  array  of  spheres.  As  we  shall  see, 
this  is  the  case  at  high  frequencies,  where  there  can  be  more  than  a  single  angle  that  the  field 
scatters.  Thus,  Eqs.  (29a),  (29b)  for  the  source  field  should  be  replaced  by  the  relations: 


Ef(r)  =  Xe-jKp"Y2R 

K 

*  {[EjJE+)(-R)e"nz  +  EUE_)(R)eT2'}ez  x  K) 
+  [Ek  M+)(~R)e  -Y2  z  +  E(™_)(R)eY2Z]K. 

-  [E(K™+)(-R)e_Y2Z  -  EKM_)(R)er2Z]— ez} 

Y  2 


where 


H(,S) (r)  =  e-jK  p-^R {-[E(KTE+)(-R)e~Y2Z  -  E(KE_)(R)eY2Z]-^K 

jcO|a0  K  ^2 

+  [E|JM+)(-R)e-Y2Z  -  E(JM_)(R)eY2Z  1— (ez  x  K) 

Y2 

+  [E£E+)(-R)e'Y2Z  +  E(JE_)(R)eY2Z  ]— ez}. 


Y  2  = 


K 


+  (72 


(57b) 


(58) 


We  have  omitted  the  subscripts  p,q  in  Kpq=kp+bpq  to  simplify  the  notation.  Here,  K  is  the  unit 
vector  of  K  and  K  is  its  magnitude.  Also,  the  amplitudes  of  the  source  field  are  given  at  z=R  and 
Z=-R,  since  our  aim  is  to  compute  the  scattering  matrix  in  relation  to  the  planes  at  z=-R  and  z=R. 
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The  source  field  above  is  in  the  plane  wave  representation.  In  the  spherical  wave  representation, 
for  each  wavevector  K,  the  source  field  can  be  writen  as  follows  (cf.  Eq.  (1119)  and  Eqs.  (34a), 
(34b)) 


E®(r)  -  2>‘JK"s'  jl[A!“S)(K)f,(a2rst)Y,S(ar!1 ) 

st  lm  (59a) 

+  A™  (K) — V  x  (f,  (02rst «  (£5r  ))], 

J<*2 

H(2St(r)  =  W(nr„) 

J»M0  „  l.  (59b) 

+  A{“S>(K)— V  x  (f,<a2ra)Y,;(nrsl ))]  , 

J°2 

where 

Affs,(K)  = - ^S-iT-e-^  (-l)1""*1 

[1(I  +  D]/2 

*{cir'[ElKTE*l(-R)  +  (-l)l‘""iE'KTE,(R)]  (60a) 

+  ^j  C{™,[e<k™*,(-R)  +  (-1)'“"*IE<k™-,(R)]}. 


a£s’(K)  =  - 


47t 


<y->  - 


[1(1  +  \)Y2  ^ 


y2R 


(-D 


l-m+1 


•{cirlErwH-D'-'Erw] 

-  Ci™[E<™"’(-R)  +  (-1)'-”*'  E'™-’(R)] } . 


(60b) 


c[^  are  given  by  Eqs.  (36a),  (36b),  where  k  is  replaced  by  K,  i.e.,  <pk  is  the  azimuthal 

angle  of  K,  and  y2  is  given  by  Eq.  (58)  Since  bpq'ast  is  an  integral  multiple  of  2n,  exp(-jK  ast)  can 
be  replaced  by  exp  (-jkp-ast)  in  Eqs.  (59a),  (59b). 

* 

For  each  wavevector  K,  multiplication  of  Eqs.  (59a),  (59b)  by  E^|n)  (rpq)  and  integration 
over  the  volume  Vpq  of  the  sphere  centered  at  apq  yields  the  following  relation 


j8g(r).E'i:l'(t„)d\xit'j‘''"A!,“s,(KK„ 


(61a) 


e_jk  p  apq  a  j^s>  (K)^, . 


(61b) 


J  Hgit(r).i®i,,)*(rpq)d3xpq 


vpq 


<?2 

j®Eo 


From  the  Ewald-Oseen  extinction  theorem  for  the  electric  field  (Eq.  (25))  it  follows  that  the 
sum  of  the  right  hand  sides  of  Eqs.  (53a)  and  (61a)  is  equal  to  zero. 


Similarly,  from  the  Ewald-Oseen  extinction  theorem  for  the  magnetic  field  (Eq.  (27)),  it  follows 
that  the  sum  of  the  right  hand  sides  of  Eqs.  (53b)  and  (61b)  is  equal  to  zero.  Since  this  is  valid 
for  all  points  in  the  lattice,  i.e.,  for  all  (p,q),  we  conclude  that  P=kp.  We  define  the  unknown 

state  vectors  B(M)(K),  B<E)(K)  with  components  b[^(K),  b|^(K),  respectively.  We  define  the 

matrices  Q(I)(kp),  Q(2)(kp)  with  components  Q^;lw(kp),Q{2).|W(kp)  and  the  vectors  A(MS)(K), 

A(ES)(K)  with  components  A[^S)(K),  A(ES)(K),  respectively.  Notice  that  B(^),B{E) are  related 

to  A,,1^,  a{E)  by  means  of  Eqs.  (46a)  and  (46b).  Finally,  we  define  the  diagonal  matrices 


r(M) 

t(M)  Pj _ £  2 

Imil'm'  _  (M) 


a 


mm'’ 


(62a) 


q(E) 

t(E)  Pj _ 

1lm;l'm'  (E) 


5,  ,6 


irumm' 


(62b) 


Then,  for  each  wavevector  K  of  the  source  field,  the  state  vectors  B(M)(K),  B(E)(K)  are  computed 
from  the  following  linear  algebraic  system: 


[i  -  T(M)Q(l)(kp)]B(M)(K)  -  T(M)Q(2)(kp)B(E)(K)  =  T(M)A(MS)(K), 
-T(E)f2(2)(kp)B<M)(K)  +  [l  -T(E)Q(l)(kp)]3(E)(K)  =  T(E)A(ES)(K), 


where  I  is  the  identity  matrix.  The  series  in  Eqs.  (54a),  (54b),  are  independent  of  the  origin  apq  in 
ap,q,  =  ap-q.  -  apq,  and,  therefore,  we  have  the  relations: 


1 


Si  p’q' 


-jkc 


U|'m'.|m(a 


p  q 


(64a) 


Q 

The  linear  algebraic  system  above  is  almost  the  same  to  that  in  the  paper  by  A.  Modinos  [5]. 
He  is  using  Bessel  and  Hankel  functions,  while  we  are  using  modified  Bessel  functions.  The 
main  difference  is  in  the  functions  U  ,m.,w  (r0),  V  lm;iw  (r0),  (Eqs.  (55a),  (55b)),  which  in  our 


(2)  rk 'v-J_Yyjkp  ap'£|' 

lm;rm'VKpJ  ~rZu 
St  p'q' 


blrm';lm(®p'q')" 


(64b) 


case  are  expressed  in  integral  form.  This  will  allow  us  to  transform  Q(l)(kp),  Q(2)(kp)  from  the 
real  domain  into  the  spectral  domain.  But  the  uniqueness  of  the  solution  of  Maxwell’s  equations 
ensures  that,  in  real  space,  the  solution  to  the  linear  algebraic  system  above  gives  the  same  fields 
to  those  of  A.  Modinos’  paper. 

Next,  we  consider  the  scattered  fields  outside  the  spheres,  i.e.,  when  |r-apq|>R  for  all  (p,q).  In 
this  case,  for  each  wavevector  K  of  the  source  field,  the  scattered  electric  field  is  equal  to 

Eg  (r)  =  2e“jkp,i,qjX  [fCWEST  VM)  +  Bj®  (IQEgrVp,)]-  («) 

pq  Im 

If  we  use  the  identity  (III  23)  in  Appendix  III,  where  we  set  o=  02,  we  obtain  the  relation 


where 


E2  k  (r)  =  ]  |  dXx^Vjx-p^2|z| 

’  2n  J  J  r T-v. 


^272 


Im 


Ml, 


(signz)— ,(px 

CT, 


1 


+  B{,^(K)- — (-jX  -  (signz) y2ez)  x  Y,), 
Jtf2 


(signz)— ,cpx 


*^e-j(kp~X)apq 

pq 


(66) 


Y2 


=  X2  +  G 


Vi 


(67) 


Here,  Y,^(w,cpx)  are  the  complex  vector  spherical  harmonics,  where  Yrm_p(Qr)  is  replaced  by 
Yrm_p(w,<Px)  *n  Eq.  (I2a)  in  appendix  I.  Also,  the  scattered  field  is  computed  outside  the 
spheres.  Therefore,  we  must  have  |z|>R. 

Eq.  (66)  can  be  rewritten  in  the  spectral  domain  by  using  the  transition  relation 

£e-j(kp-x).apq  ^^£§(2>(X_kp-bp,q,),  (68) 

pq  p'.q' 


where  Sc  is  the  area  of  the  unit  cell  in  real  space  and  bpq  has  already  been  defined  (Eqs.  (48a), 
(48b),  (49)).  The  scattered  electric  field  becomes  equal  to 


Eg(r)  =  jvErV^K'P',2l,T 

Sc  K'  a2Y2  1m 


B|m)(K)Y1^f((signz)— ,(pK.) 


+  B{E)(K)-t^— (-jK'  -  (signz)y'2ez)x  Yj™ ((signz)-^-,q>K>  ) 


J<*2 


where 


72 


K'  +o. 


1/ 

/2 


(69) 


(70) 


We  have  omitted  the  subscripts  p',q'  in  Kp.q.  =kp  +bp-q.  to  simplify  the  notation.  A  straight¬ 
forward  computation  leads  to  the  following  two  relations 


where 


[1(1  +  Yjfi  Y,;((signz)2i<pK)  =  (signz)'-m*l[c(™>K  +  C£E»(ez  x  K)] 
a2 

+  (signz)'~mD{™)ez 


(-jKK-  (signz)y  2ez  )x  [l(l  +  #  Y/" 


(signz)— ,(pK 

a2  J 


-(signz)'-m+1jKCSE)ez 


p(TM)  _„my 
Ulm  ~a' 


l,m+l 


Y2_ 

VCT2 


■<Pk 


,-J‘PK 


m-1 


y_2_ 

Va2 


■<Pk 


J<PK 


r(TE)  --i 
Mm  “  J 


- 

r  \ 

p-jl>K  _RmV 

e  Pi  Yl.m-1 

—  ><Pk 

ej(PK 

\a2  ) 

lCT2  J 

. 

or 


=  mY, 


lm 


y2 


Va2 


-<Pk 


Y  2 


(71a) 


(71b) 


(72a) 


(72b) 


(72c) 


Here,  K  is  the  unit  vector  and  K  the  magnitude  of  K.  Substitution  of  Eqs.  (71a),  (71b)  into  the 
scattered  electric  field,  yields  the  following  relation  in  the  plane  wave  representation,  when 
|z|>R: 
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where 


E(2  kM  -  Ee-iR' "  ®(z  -  R)e-'i<l-R»  x  K')+  E^^R)!  K'  - J^e, 

K’  [  L  V  Yl 

+  ©(-z-R)eT2(z+R|  E£,(-R)(i,xK,)tEg;->(-R/k'  +  £t,|  j 


,  (73) 


Eg'W^-^Y  1  -j7[-!-rcglBiy>(K)-4cr>(K)Bgl(K)1  (74a) 


3c  lm 


[l(l  +  l)]/2La2Y2 


'[l(l  +  l)P  L®2Y2 


-2^  C<™'b!“»(K)  +  -L  Cj^>BjR>(K)l ,  (74b) 


Effg’f-R)  -  ^e-'2»  ■£  y  1)'"mii7  *  [ -Acff  >Bj»>(K)  +  4rCi™>Bim (K)l ,  (74c) 

C  [l(l  +  l)p2  LCF2Y2 


EE-»(-R)  =  ^e-'2RX^ 


Sc  ta[l(l  +  l)]>2  La2Y2 


*  -2-7c!mM,B|m’(K) — VCi™>B!m  (K)  .  (74d) 


The  scattered  magnetic  field  is  computed  from  the  curl  of  the  electric  field  and  is  equal  to 


Hjk (r)  =  t^2- Z e_jK'P  ®(z - R)e 
J«Eo  R'  [ 


-Y'2(z-R)  a2  p(TM+) 


fE^c+'(R)(ezxK') 


-E(kt,esc+)(R)  — K'  ez  +©(-z-R)e^(7-+R)  -^-E^C-RXe,  x K') ,  (75) 

(a2  cr2  J  i  Y2 


+  ES£>(-R)f22.K'  +  i51«1 


a2  g  2 


At  this  point,  we  introduce  the  scattering  matrix  that  relates  the  fields  at  z=-R  and  z=R  by 
means  of  the  relation 


'e£-’(-R) 

e<™;»(-r) 


E<kte*»(-R)' 

E(K™+)(-R) 


ee*(r)  =Sk"  Er’(R) 

eS’(R)  e'™-'(R) 
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The  components  of  the  right-hand  vector  are  the  amplitudes  of  the  source  electric  field  (of  Eq. 
(57a))  while  the  components  of  the  left-hand  vector  are  the  amplitudes  of  the  scattered  electric 
field  (cf.  Eqs.  (73),  (74a)-(74d)).  The  matrix  elements  of  SKK  are  obtained  by  solving  the 
algebraic  system  in  Eqs.  (63a),  (63b)  for  each  of  the  four  individual  cases  where  one  of  the 
components  of  the  right-hand  vector  in  Eq.  (76)  is  set  equal  to  one,  and  the  rest  are  zero.  Then, 

in  each  case,  we  use  the  solutions  for  B(M^(K),  B^E)(K)  to  compute  the  scattered  amplitudes 
from  Eqs.  (74a)-(74d)  and  identify  them  as  the  elements  of  SKK . 


If  we  divide  SK1C  into  submatrices,  i.e., 


S 


K'K 


o(2,l) 

^K'K 


(77) 


then  S<">  S<'£>  give  the  (TE)  and  (TM)  reflected  and  transmitted  amplitudes  of  the  electric  field 
at  z=-R  and  z=+R,  respectively,  with  wavevector  K' ,  if  the  source  field  is  either  a  forward  (TE) 
or  (TM)  wave  at  z=-R  with  wavevector  K.  Similarly,  give  the  reflected  and 

transmitted  amplitudes  at  z=+R  and  z=-R,  respectively,  if  the  source  field  is  a  backward  wave  at 
z=+R.  We  see  then  that  the  scattering  matrix  SK  K  contains  all  the  information  about  the 
scattering  of  a  plane  wave  by  the  array  of  dielectric  spheres. 

The  total  field  at  z=-R  and  z=R  is  the  sum  of  the  source  field  and  the  scattered  field.  If  we 
define,  for  the  total  field,  the  two-component  vectors 


E(k+)(r)  = 


E(kE+)(R) 

E(™+)(R) 


(78a) 


E(k)(R) 


E(kte_)(R) 

E(k™-)(R) 


(78b) 


and  similar  vectors  E<^)(-R),  E^^-R)  at  z=-R,  then  it  follows  from  Eqs.  (57a),  (76),  (77)  that 

Efc’(-R)  =  SilliE>K+»(-R)  +  (e-2«RSK.KIM  +  SK>)e«k-»(R)  .  (79a) 

E£>(R)  =  (e-2,2R8KKI2.2  +  S<“  jEir>(-R)  +  s£-i>Efc>(R) ,  (79b) 

where  I2x2  is  the  2x2  identity  matrix  and  5KK  =1  if  K'  =  K  and  zero  otherwise.  Here,  we 
added  the  amplitudes  of  the  source  field  to  the  transmitted  components  of  the  scattered  field.  In 
matrix  form,  Eqs.  (79a),  (79b)  can  be  written  as  follows: 
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E(k-)(-r)1  =  r  S<^(-R,R)  Sgg(-R,R)I: Et>(-R) 
E(k+>(R)  I  [sgi)(-R,R)  S^(-R,R)JL  E(k->(R) 


(80a) 


where 


o(U)  /T»  D\_e(M) 

i»K,K(  k,k;-&k,k, 


(80b) 


S&2(-R,R)  =  e-2^R8K.KI2)<2  +  S<‘-.2> , 


(80c) 


(-R’R)  =  e  2^2R^K'K^2x2  +Sk:k  > 


(80d) 


q(2,2)/_p  n\_c(2,2) 


(80d) 


In  the  relation  above,  both  the  right-hand  and  left-hand  vectors  refer  to  the  amplitude  of  the  total 
electric  field  at  z=— R  and  z=R. 

Given  the  total  scattering  matrix  SK-K(-R,R),  it  is  easy  to  obtain  the  total  scattering  matrix 

SK'K  (-zbz2)  that  refers  to  the  amplitudes  at  z=-zi  and  z=Z2,  where  zi>R,  Z2>R,  using  the 
propagation  matrices.  From  Eqs.  (57a),  (73),  we  see  that 


Ek^(z2)  -  Er (z2  < - zi)Ek^(-zi)  ’ 


(81a) 


E(k)(z2)  =  PkI(z2  ♦--zOE^C-zi), 


(81b) 


where 


Pr(z2  * - zl)~e  Y2<'Z2+Z'^2x2’ 

Pk*(z2  < — Z])  =  e1'2(Z2+Z|)l2x2- 

Substitution  of  Eqs.  (81a),  (81b)  into  Eq.  (80a)  leads  to  the  relation 


(82a) 


(82b) 


ek;>(-zi)1  c  /  K+)(-zi) 

Er'  (z2)  J  LEk(z2) 


(83a) 


where 


^K'k(  zhz2)- 


^K'k(  zbz2)  (  Z„Z2) 

§K'K  (-ZI’Z2)  S£2)(-Z,  ,Z2)J 


(83b) 
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and 


• 

(-Z,  >  z2)  =  e'ty2+T2Kzi  -R)Sj^(-R,  R) , 

(84a) 

sK,(-z„z2)  =  e-rt(,'-R>-,2<,2-R,sK»(-R,R), 

(84b) 

• 

Sr'k  (-zi>z2) =  e  ^,Z2_r)”t2<Z|~R)Sk-k  (-R,R), 

(84c) 

S&2,(— zi.z2)  =e  «-^y2"z2-R»sg'(-R,R) . 

(84d) 

• 

The  transfer  matrix  T(z2  < — z,)  can  be  easily  obtained  from  Eqs.  (83a),  (83b)  by  constructing 
the  state  vector  that  refers  to  all  the  wavevectors  K  under  consideration.  We  group  together  all 
the  backward  waves  and  all  the  forward  waves  and  define  the  four  state  vectors 

(85) 


where  a  =  (+)  or  (-) ,  i.e.,  forward  or  backward  wave  (of  Eq.  (78a),  (78b))  and  z  =  -z\  or  z2. 
Then,  it  follows  from  Eqs.  (83a),  (83b)  that 


E<->(-Zl)  =  S^E^-z,)  +  S(I’2)E<2->(z2)  , 


(86a) 


E(2+)  (z2  )  =  S(2’1)E(2+)(-z1  )  +  S(2’2)E^_)  (z2  ) ,  (86b) 

where  the  matrix  elements  of  S^1’^  are  equal  to  =  S£’^Kn  (-z,,z2),  m=l,2,...,N, 
n=l,2,...,N.  The  transfer  matrix  is  computed  from  the  algebraic  system  of  Eqs.  (86a),  (86b),  and 
expressed  by  the  relation 

E2(+)(z2) 

_E(2_)(z2). 


=  T(2 


'2  < - Zl 


(-Zl) 


(87a) 


where 


T(z2  < — Z|) 


g(2,l)  _g(2,2)g(l,2r1g(l,l)  g(2,2)g(l,2)  1 

_g(l,2)-|g(l,l)  S(Ur’ 


(87b) 


•  The  computation  of  the  transfer  matrix  allows  us  to  superimpose  many  layers  of  dielectric 

spheres  at  the  same  or  different  distances  between  them  and  compute  the  transfer  matrix,  as  well 
as,  the  scattering  matrix  of  the  total  set  of  layers. 


23 


The  wavevectors  Kpq  that  should  be  included  in  the  state  vectors  above  (cf.  Eq.  (85))  are  for 
all  the  propagating  modes.  In  a  non-absorptive  host  medium,  where  the  real  part  of  ct2  tends  to 
zero,  a  pair  of  integers  (p,q)  for  which  Y2(Kpq)  (cf.  Eqs.  (58),  (70))  is  purely  imaginary  (in  the 

limit  where  Re(a2)  -*  0),  defines  a  propagating  mode  in  medium  2  (cf.  Eqs.  (84a)-(84d)).  In 
particular,  the  mode  for  the  source  field  for  which  p=q=0,  is  the  principal  propagating  mode. 
The  rest  of  the  modes  are  evanescent.  In  addition  to  the  fact  that  the  amplitudes  associated  with 
them  decay  exponentially  with  distance  (cf.  Eqs.  (84a)-(84d)),  their  radiated  power  is  equal  to 
zero,  i.e.,  they  do  not  contribute  to  the  reflection  and  transmission  coefficients.  Therefore,  there 
is  no  need  for  them  to  be  included  in  the  state  vectors.  On  the  other  hand,  in  an  absorptive  host 
medium  where  cr2 ,  has  a  positive  real  part,  we  define  the  propagating  modes  as  the  pairs  (p,q) 
for  which  y2(Kpq)  satisfies  the  criterion  |lm(y2)/y2|  >  e .  Here,  s  is  a  small  positive  number, 

e.g.,  £  =  0.001 ,  that  should  be  adjusted  so  that  the  rest  of  the  modes  have  a  negligible 
contribution  in  the  radiated  power  by  the  array  of  dielectric  spheres. 

IV.  n|»lw(r0),Q£>lw(r0)  in  the  real  and  spectral  domain 


The  functions  O[^.|,m.(r0),  Q{^rm.,(r0)  are  given  by  Eqs.  (64a),  (64b)  and  the  functions 

U|m;iw(ro)>  V,m;i'm’(ro)  are  given  by  Eqs.  (55a),  (55b).  In  Appendices  IV  and  V  the  scalar  and 
vector  addition  theorems  are  derived.  Thus,  from  the  vector  addition  theorem,  i.e,  Eq.  (VI),  we 
obtain  the  relation 


F(M,out) 

&2M 


(r  -  r0)  =2[u,ra;,'m.(r,)E«».;»(r)  +  V,„,,m.(r0)Eg;>(r)J, 


I'm' 


(88) 


where  Ulm.rm.(r0),  V)m;iw(r0)  are  given  by  Eqs.  (V2a),  (V2b),  and  a  =  cr2  in  these  equations. 
Similarly,  from  Eq.  (V3),  are  obtained  the  relation 

Eyr‘Hr-r0)  =  Xl°l».Jv(ro)E!.j™”)(r)  +  Vlm;|.m'<>-o)ES,<r)J.  (89) 

lm 

Substitution  of  Eq.  (88)  into  Eq.  (55a)  and  Eq.  (89)  into  Eq.  (55b),  together  with  the 
orthonormality  of  the  vector  spherical  harmonics  and  the  definition  of  C,\  (cf.  Eq.  (56))  lead  to 
the  relations 


U|m;l'm'(ro)  —  ClUim;l'm'(ro)  »  (90a) 

(ro)  =  ^|V|m;|'m'(ro).  (Mb) 


Therefore,  in  the  real  domain,  ^im.|'rn'(ro),(iim  i'm'(ro)  (cf-  Eqs.  (64a),  (64b))  are  given  by  the 
following  relations: 


nS'Jm(kp) 

\ 


[l(l+l)l'(l’+l)]X 

.2ajnP?'+1Zi;£;,(kp)  +  2a 


*  I9«mft!1,'+I7!'’m'+Vb  1 4. 9n»j!1'_1(3jnzj,I^7,(kp)  +  mm'Zj^1  (kp)] 


(91a) 
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(91b) 


where 


4Jilm(kp)  = 


_ 2T+1 _ (_1) 

[l(l+l)l'(l'+l)]/7 


m  /  471 

V  3 


»[V2arB,._1.m,tl(l,-l;rm')Z!;;l;7*l(kp) 
+  mB,-.,  m  (l,0;rm')Zj  m  (kp)] 


Z!mm'(kp)  =  S'e' 


and  Glm;l'm'(apq)  is  Siven  bY  Ecl-  (IV  2)>  where  a  =  °2  • 

These  expressions  are  identical  to  those  in  the  paper  by  A.  Madinos  [5]  except  for  the  function 

Gim  i'm'(apq)  since  be  is  using  Hankel  functions  while  we  are  using  modified  Bessel  functions. 

We  see  then  that  the  integral  equation  approach  provides  the  same  solution  as  the  approach  of 
applying  the  boundary  conditions  at  the  surface  of  the  dielectric  spheres. 

The  computation  of  Q{^;1,m,(kp),  Qgrm,(kp)  in  the  spectral  domain  requires  the  Fourier 

transforms  of  the  vector  functions  in  Eqs.  (55a),  (55b),  which  are  derived  in  Appendix  VI.  In 
terms  of  the  Fourier  transforms,  i.e.,  Eqs.  (VI  8a)  -  (VI  8c),  Eqs.  (55a),  (55b)  become  equal  to 


flm;l'm,(ro)  (2jf)3 


.3 

( 2n)J 


jE(2EJt,(X)-E(2^n)*(X)eJ"r°d3X, 


(93a) 


(93b) 


the  integration  extending  over  the  whole  space  in  the  spectral  domain.  Notice  that,  when  r0  =  0 

in  Eqs.  (55a),  (55b)  then  the  functions  E^i’^^r) ,  where  a  =  M,E,  are  zero  inside  V0  and  the 

functions  U  im,rm'  (0),Vlm;lW(0)  are  zero.  Therefore,  when  the  integral  form  of  U  and  V  is 
used,  then  the  prime  in  Eqs.  (64a),  (64b)  can  be  removed,  i.e.,  the  term  with  p=q=0  can  be 
included  in  the  summation.  Thus,  substitution  of  Eqs.  (93a),  (93b)  into  Eqs.  (64a),  (64b)  and 
application  of  the  transition  relation  in  Eq.  (68),  leads  to  the  relations 


n!''f*(kp)=zriZi  KS’ 


(94a) 


nS-j»(kp)  =  if  srEi  KS’V.x.)  •  Egr ’(K.x.xu, 


(94b) 


where  Kpq  =  kp  +  bpq ,  and  the  summation  is  over  all  (p,q)  points  in  the  reciprocal  lattice.  The 
subscripts  (p,q)  have  been  omitted  in  Eqs.  (94a),  (94b)  for  simplicity  of  notation. 
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Next,  we  compute  the  right  hand  side  of  Eq.  (94a).  From  the  definition  of  the  vector 
spherical  harmonics  (cf.  Eqs.  (I la),  (V9)),  we  obtain  the  relation 

Y^(QX)-Y^(QX)  =  - - 

[l(l+l)l  (1  +l)]/2 

*  [2of  Y,:  m,+l  (Q,  )Y,m+i  (QJ  +  2P?'prY, (nOY,.™-,  (□,.).  (95) 

+  mm'Yi:m.(Qi)Ylm(t!).)] 

where  X  =  (kp,Xz)  and  A,p=Kpq.  Substitution  of  Eqs.  (IV  8a),  (VI  8b),  (95)  into  Eq.  (94a) 
leads  to  the  relation 

q(|)  tk  )  = _ ! - 

lm'lnt  p)  (KltDTO'-tDlK  (96) 

*  [2arpr*'z!:™^'(kp)  +  2aj’,-|przi:c:i'(kp)  +  mm'Z|^(kp)]  , 

where 


yl’m'/K  \  _  87lR^ 

^Im  lKpl-  scQ' 


•^-(^02)  G|(X..ct2) 
^1+Y22  xz+Y2 


•C2=K2+CT2  (98) 

The  functions  F1(X.,CT2)5Gl(^’a2)  are  §iven  bY  Ecls-  (VI  3  b),  (VI  4b),  respectively.  Also, 
A.  =  (x2  +  X\ Y2 ,  where  Zp  =  Kpq  and Kpq  =  kp  +  bpq ,  bpq  being  a  vector  in  the  reciprocal 

~  |  >  t 

lattice.  Notice  that  Eqs.  (91a),  (96)  are  identical  except  for  Zlm  (kp)  being  replaced  by 
zj  ™  (kp ) .  It  is  rather  obvious  then  that  we  should  have  the  identity 

z!mm'<kp)  =  Z!mm'<kP>-  <99) 

In  Appendix  VII  it  is  shown  that  this  is  the  case.  Eqs.  (96),  (97)  provide  Qj!^,.im(kp)  in  the 
spectral  domain  and  Eq.  (91b)  provides  fl|^,.Im(kp)  in  the  spectral  domain  if  zj™  (kp)  is 
replaced  by  z|^'(kp). 

The  integral  in  Eq.  (97)  can  be  computed  analytically  only  when  K=0,  i.e.,  for  normal 
incidence  of  the  source  wave  and  at  the  origin  of  the  reciprocal  lattice.  In  the  general  case,  the 
integral  can  be  expressed  as  a  series  of  terms  involving  modified  Bessel  functions.  The  explicit 

i  r  t 

expression  of  the  series  will  be  given  elsewhere.  Therefore,  the  computation  of  Zj^  (kp)  in 
Eq.  (97)  involves  a  three-dimensional  summation,  as  compared  to  a  two-dimensional  summation 
for  the  computation  of  (kp)  in  Eq.  (92).  On  the  other  hand,  for  large  values  of  (p,q),  the 

i  i-3 

sum  in  Eq.  (97)  is  of  order  Kpq  (after  the  integral  has  been  expressed  as  a  series),  while  the 
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sum  in  Eq.  (92)  is  of  order 


lpq 


-1 


Therefore,  in  a  lossless  medium,  the  sum  in 


7l'm' 

Alm 


(kp) 


should  converge  much  faster  than  that  in  zj  ™  (kp),  notwithstanding  the  third  sum  due  to  the 

integral  in  Eq.  (97).  Due  to  the  oscillatory  behavior  of  the  terms  in  Eq.  (92),  the  maximum 
values  of  (p,q)  must  be  at  least  equal  to  104  in  order  to  obtain  the  plot  in  Fig.  3  in  the  paper  by  N. 
Stefanou  [6].  Therefore,  4*108  terms  are  required  in  the  sum  of  Eq.  (92)  for  convergence.  On 
the  other  hand,  some  preliminary  results  using  Eq.  (97)  indicate  that  the  maximum  values  of 
(p,q)  must  be  equal  to  50  and  the  number  of  terms  in  the  sum  for  the  integral  should  be  200. 
Therefore,  there  are  5*105  terms  to  be  summed  in  Eq.  (97).  Thus,  there  is  an  advantage  to 
perform  the  computations  in  the  spectral  domain  rather  than  in  the  real  domain  for  lossless  host 

~  j  f  t 

medium.  On  the  other  hand,  in  a  lossy  host  medium,  for  sufficiently  high  absorption  Z,m  (kp) 


will  converge  faster  than  zj^1  (kp)  due  to  the  exponentially  decaying  terms  in  the  modified 
Bessel  functions  (cf.  Eqs.  (92),  (IV9),  (I7b)). 


V.  Multilayer  and  Three-Dimensional  Array  of  Spheres 

When  there  are  more  than  one  layer  of  dielectric  spheres,  each  layer  is  obtained  from  the 
previous  one  by  shifting  all  the  spheres  of  the  latter  layer  by  the  vector  a3 .  At  the  zeroth  layer, 
which  lies  at  z=0,  if  the  spheres  are  located  at  apq,  then  at  the  nth  layer  the  spheres  are  located  at 
apq  +  na3.  The  translation  of  the  spheres  from  apq,  at  z=0  to  apq+napq  is  equivalent  to  moving  the 
coordinates  of  the  nth  layer  from  r  to  r-na3.  Therefore,  for  the  nth  layer,  instead  of  Eq.  (73),  the 

1  th 

scattering  amplitudes  at  z  =  ±  A a3z ,  where  z  is  measured  from  the  n  layer,  become: 


a3z) 


=  e 


jK’na3-Y2(Fa3z-R) 


E&cTO , 


(100a) 

(100b) 


Also,  for  the  n,h  layer,  instead  of  Eq.  (57a),  the  amplitudes  for  the  sourse  field  at  z  =  ±^a3z  are: 


E<K,(n;-i‘>,I)  =  eJK”3*’2^*3^RW’(-R). 

EK,(n;ia3z)  =  eJ 


(101a) 

(101b) 


On  the  left  hand  side  of  the  relations  above,  the  first  index  indicates  the  layer  and  the  second 
index  is  the  value  of  z  (measured  from  the  n,h  layer)  where  the  amplitudes  are  computed.  By 
adding  the  amplitudes  of  the  source  field  to  the  transmitted  components  of  the  scattered  field 

from  z  =  --^-a3z  to  z  =  ^a3z,  instead  of  Eqs.  (80a)-(80d)  we  obtain  the  following  relation  for 

the  scattering  matrix  of  the  nth  layer: 


2  a3z)  _  ej(K'-K)  na3 
pW/tvir,  'k 


^K' (n;  2a3z)  J 

‘sK(-ia3zla3z)  Sgi,(-ia3Zla3Z)lE«(n;-la3z) 
SjSHaa. ,±a3z)  S^(-ia3z  ,{a3z)j[  E^n;^.) 


(102) 


where  sj^(-zl’z2)  is  given  by  Eqs.  (84a)-(84d).  If  we  define  vectors  similar  to  those  of  Eq. 
(85)  for  the  nth  layer  and  also  define  the  diagonal  matrix 


G(A)  = 


HK2-a 


(103) 


'•e-jKNA 


then,  instead  of  Eqs.  (86a),  (86b),  we  obtain  the  relation 

rG(na3)E(->(n;-ia3z)1  [§0.0  §^2)  f  G(na3)E«(n;-la3z) 


G(na3)E(+)(n;Ia3z) 


:(2,l)  c(2,2) 


G(na3)E«(n;ia3z) 


(104) 


The  transfer  matrix  is  computed  from  the  relation  above  for  the  nth  layer,  from  the  plane  at 
z  =  “a3z  to  the  plane  at  z  =  ^a3z ,  namely. 


where 


E(n  +  1)  =  G(-na3)T(^a3z  4-  -|a3z)G(na3)E(n), 


'E<+>(n;-ia3z) 

E(n)  =  /  3  f 

EG)(n;-la3z) 


(105) 


(106a) 


E‘+>(n;ia3z) 
E(n  +  1)  =  ,  f 

E  }(n;^a3z) 


G(A)  = 


G(A)  0 
0  G(A) 


(106b) 


(107) 


and  T(z2  < — Z| )  is  given  by  Eq.  (87b). 

For  n  layers,  it  follows  from  Eq.  (105)  that 


E(n)  =  G(-(n  -i)a3)TnG(-la3)E(0) , 


where 


(108) 


T  =  G(2-a3)T(-f-a3z  4-  -±a3z)G(^a3) . 


(109) 
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Eq.  (108)  can  also  be  written  as  follows: 

E(n)  =  TnE(0) ,  (110) 

where 

E(n)  =  G((n~i)a3)E(n).  (Ill) 

Then,  for  two  consecutive  layers  we  have  the  relation 

E(n  + 1)  =  f  E(n) .  (112) 

For  a  three-dimensional  lattice,  i.e.,  when  there  is  an  infinite  number  of  layers,  application  of 
the  Bloch  theorem,  namely, 

E(n  +  l)  =  ejkza3zE(n),  (113) 

leads  to  the  homogeneous  algebraic  system 

[f-ejkza3zIjE(n)  =  0.  (114) 

The  dispersion  relation  of  the  three-dimensional  lattice,  i.e.,  kz  as  a  function  of  (kx,ky,oo)  is 
determined  by  the  determinant  in  Eq.  (114)  being  set  equal  to  zero.  Following  this  procedure, 
we  can  compute  the  propagating  modes,  as  well  as  the  band  gaps  in  the  infinite  lattice. 

Another  approach  to  determine  the  disperson  relation  of  the  infinite  lattice  is  as  follows:  Let 
the  unit  cell  of  the  three  dimensional  lattice  be  determined  by  the  translation  vectors  ai,  a2,  a3,  so 
that  the  lattice  points  are  given  by  apqr=pai+qaz+ra3.  The  fundamental  vectors  bj,  b2,  b3  of  the 
reciprocal  lattice  are  given  by  the  usual  expressions  in  term  of  aj,  a2,  a3  and  they  satisfy  the 
relations  a;  bj  =  2^ .  Notice  that,  if  a3  -ez  *  0 ,  then  bi,  b2  are  not  parallel  to  B|,  B2  in  Eqs. 

(48a),  (48b).  The  points  in  the  reciprocal  lattice  are  given  by  bpqr  =  pb,  +  qb2  +  rb3 .  Following 

the  same  steps  in  Section  III,  we  conclude  that  in  Eqs.  (64a),  (64b)  become  functions 

of  the  three-dimensional  vector  k,  the  sums  are  over  (p,q,r)  and  kp,  apq  are  replaced  by  k,  apqr, 
respectively.  The  dispersion  relation  of  the  lattice  is  determined  by  setting  the  right-hand  side  in 
Eqs.  (63a),  (63b)  equal  to  zero,  and  computing  the  values  of  kz=f(kx,  ky,  co)  for  which  the 
determinant  of  the  left-hand  side  is  equal  to  zero.  The  functions  Q(1)  (k),  Q(2)  (k)  are  still  given 
by  Eqs.  (91a),  (91b),  (92),  where  the  sums  are  over  (p,q,r)  and  kp,  apq  are  replaced  by  k,  apqr, 
respectively.  Thus,  Eq.  (92)  becomes 

zLm'(k)  =  Ze-i,l-'''Glm;1.m.(apqr).  (115) 

pqr 


where  G|m;,.m.(apqr)  is  given  by  Eq.  (IV  2)  and  a  =  a2  .  Also,  Eq.  (99)  is  still  valid,  except  that, 
in  the  three-dimensional  case,  Eq.  (97)  becomes  equal  to 


(116) 
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Here,  Vc  is  the  volume  of  the  unit  cell  and  Kpqr  =  k  +  b 


pqr ' 


that  (k)  is  of  order 


‘pqr 


,  while  Zj  ™  (k)  is  of  order 


In  a  lossless  host  medium,  we  see 

1-4 


K 


pqr 


Therefore,  the  sum  in  Eq. 


(1 16)  converges  much  faster  than  that  in  Eq.  (115).  On  the  other  hand,  in  a  lossy  host  medium. 


for  sufficiently  high  absorption,  Zj^'(k)  converges  faster  than  Zj™  (k)  due  to  the  exponentially 
decaying  terms  in  the  modified  Bessel  functions  (cf.  Eqs.  (115),  (IV  9),  (I  7b)).  It  is  only 
through  actual  computation  that  can  be  determined  whether  Eq.  (115)  or  Eq.  (116)  should  be 
used  in  the  presence  of  a  lossy  host  medium. 


VI.  The  Array  of  Spheres  Embedded  in  a  Slab 


The  slab  is  medium  2  with  complex  wave  number  c2  (cf.  Eq.  (1))  and  lies  between  the  planes 
at  z=-zi  and  z=Z2  (zi,  Z2  are  positive).  The  medium  1  for  z>Z2  has  wavenumber  01  and  the 
medium  4  for  z<-zi  has  wavenumber  04.  Here,  we  consider  the  general  case,  while,  in  most 
cases,  media  1  and  4  are  the  vacuum.  The  two-dimensional  array  of  spheres  lies  at  z=0,  and  its 
scattering  matrix  is  given  by  Eqs.  (86a),  (86b),  namely, 


E^(-z,) 
.  E(2+)(Z2) 


(117) 


For  a  multilayer  array,  the  scattering  matrix  can  be  obtained  from  Eq.  (108)  and,  therefore,  Eq. 
(117)  covers  this  case  too. 


The  vectors  E(2±)(-z1),E(2±,(z2)  are  given  by  Eq.  (85),  i.e.,  they  refer  to  all  the  wavevectors  K 
relevant  to  the  computation  (cf.  the  paragraph  at  the  end  of  Section  III). 

Application  of  the  boundary  conditions  for  the  electromagnetic  field  at  z=-zi,  z=Z2  for  each 
wavevector  K  yields  the  relations 


• 

t24E(2+)(_z,)  =  E(4+)(-z,)  +  r24E^(-z,), 

(118a) 

t24E(2_)(-z,)  =  r24E4+)  (-z, )  +  e£->(-z,). 

(118b) 

t12Ej+)(z2)  =  E2+)(z2)  +  r12E^-)(z2), 

(118c) 

• 

t12E|-)(z2)  =  r12E<+)(z2)  +  E(2-)(z2). 

(11 8d) 

Here,  ry,  ty  are  diagonal  2Nx2N  matrices  with  diagonal  elements 

rfE’(KN),  fM»(KN)  and  tf»(K,),  t<™>(K,),  ....  tl/E»(KN), 
where 

rfE)(K,),  rij™)(K1),  ..., 
t[jTM)(KN) ,  respectively, 

• 

rij  (K)-Yi+Yj’ 

(119a) 

(119b) 

« 

rr«0=frt, 

<Vj+<TjYi 

(119c) 
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t!™'(K)  = 


2„. 


2afYj 


2  2 
CTi  yj+CTj Yi 


Yi  = 


K2+a? 


(1 19d) 
(11 9e) 


Also,  the  subscript  i  in  E-±)(z),  denotes  the  medium  where  the  amplitudes  are  computed. 
•  Rearranging  Eqs.  ( 1 1 8a)  -  (1 1 8d),  we  obtain  the  following  relations 


~E(2+)(-z,)~ 
_  E2-)(z2)  _ 

=  A 

E^-z,)' 
_  E|_)(z2)  _ 

-B 

'e2-)(-Z[)" 
.  E2+)(z2)  _ 

E^(-z,)' 
_  E{+)(z2)  _ 

=  B 

"Ei+)(-z,)" 

+  D 

^(-z,)' 
_  E2+)(z2)  _ 

where 


(120a) 


(120b) 


A  = 


I  -  r24  0 


M2. 


B  = 


r24  0 


M2. 


D  = 


l24 

0  I  -r. 
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(121a) 

(121b) 

(121c) 


Here  1,0  are  the  2Nx2N  unit  and  zero  matrices,  respectively.  The  solution  of  the  algebraic 
system  of  Eqs.  (117),  (120a),  (120b)  yields  the  relation 


E{+)(z2) 


(122) 


where 


s  =  b  +  ds[i  +  bs]~'a. 


(123) 


S  is  the  total  scattering  matrix  for  the  total  amplitudes  computed  just  outside  the  slab  and, 
therefore,  it  provides  the  reflection  and  transmission  coefficients  of  the  array  outside  the  slab. 


VII.  Reflection  and  Transmission  Coefficients 

From  the  previous  section  we  see  that,  if  the  source  and  scattered  field  have  wavevectors  kp 
and  K' ,  respectively,  then  their  amplitudes  in  media  1  and  4  are  related  by  the  elements  of  the 
total  scattering  matrix  in  Eq.  (122),  namely, 
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'eS^c-zi)! 

E^-z,)  -  ,  E<™+)(-z,) 

4K  1  =  S(K  ,k  )  p  , 

E^+)(z2)  P  E[™->(z2) 


E|r'(z2) 

e!™+)(Z2) 


(124) 


ec-w 


The  power  flux  in  medium  i,  at  z  and  through  the  area  Sc  of  the  unit  cell  is  equal  to 

Pi  (z)  =  Re  J|e* (p, z)  x  H;  (p, z) j •  ezdxdy . 


(125) 


In  the  plane  wave  representation  (cf.  Eqs.  (57a),  (57b),  (73),  (75)),  the  power  flux  for  the  forward 
(+)  and  backward  (-)  waves  becomes  equal  to 


P,«»(z)  =  Sc2:  Refe)E!™(z)2  +  Re 


jconOYi(M  iK  v  ' 


(126) 


where 


fi(K')  =  [K 


'2  +  (J?  I  2. 


(127) 


Thus,  for  a  backward  source  field  in  medium  1,  the  power  flux  for  the  (TE-)  and  (TM-)  source 
waves  at  z=Z2,  are  equal  to 


PS!-’(z2)  =  ScRe(W]  |EjJ;»(z2)  , 


P.(™-)(z2)  =  ScRe 


jw^oYURp 


(128a) 


(128b) 


Similar  expressions  are  obtained  for  the  source  field  in  medium  4. 

The  reflection  and  transmission  coefficients  follow  from  Eqs.  (124),  (126),  (128a),  (128b). 
For  example,  if  the  source  field  is  a  (TE-)  wave  with  wavevector  kp  in  medium  1,  then  the  total 
reflection  coefficient  at  z=Z2  is  equal  to 


/e,(kD,TE)  =  l 


Re(-jyi) 


(129) 


•£  Refc-jy,)(8K..tp-S33(K',kp)](81[,^+S^(K'.kp))]+Rj^-  |s43(K',kp)f 
k'L  Vjyi  J 

where  y,  =  [k2  +ct2  P  ,  y'|  =  [k'2  +  c2  ]  2 .  The  total  transmission  coefficient  at  z=-zi  is  equal  to 


-^4(kp*TE)  Re(-jyj ) 


yZ  Re(-jY'4)|S,3(K',kp)  +Re  -f-  S23(K',kp) 


(130) 
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On  the  other  hand,  if  the  source  field  is  a  (TM-)  wave  in  medium  1 ,  then  the  total  reflection 
coefficient  at  z=z2  is  equal  to 


/f,(k  ,TM)  =  1-- 


Re 


(  „2' 
Zl 
in 


K' 


Re(-jr;)S,4(K'.kp) 


+  Re 


-s„(K',kp)XsK..Up +s;4(K',kp)) 


and  the  total  transmission  coefficient  at  z=-zi  is  equal  to 


(131) 


_74(kp,TM)  = 


Re! 


JY1 

V  J 


K 


Re(-jY'4)S14(K',kp) 


f  .2  \ 


+  Re 


vJY4y 


S24(K',kp) 


(132) 


where  y4  =  K'2  +  a4 J2 .  Similar  expressions  can  be  obtained  if  the  source  field  is  a  forward  (+) 
field  in  medium  4.  Notice,  that  each  individual  term  in  the  sums  in  Eqs.  (129)-(132)  provides 
information  as  to  how  much  each  individual  wavevector  K  contributes  to  the  total  reflection  and 
transmission  coefficients.  In  the  case  when  the  host  medium  or  the  spheres  between  the  planes  at 
z=-zi  and  z=z2  are  lossy,  the  absorption  coefficient  A,  is  equal  to 

A  =  \-R-J ,  (133) 


where  R.  ./refer  to  each  case  mentioned  above  for  the  reflection  and  transmission  coefficients. 

VIII.  Summary  and  Conclusions 

A  semi-analytic  integral  equation  approach  was  presented  in  this  paper  for  the  computation 
of  the  reflection  and  transmission  coefficients  by  a  single  or  multi-layer  array  of  dielectric 
spheres  embedded  in  a  host  medium.  It  was  shown  that  this  approach  is  equivalent  to  that  of 
applying  the  boundary  conditions  at  the  surface  of  the  spheres.  The  integral  equation  formalism 
is  simple  in  its  application  only  for  special  shapes  of  scattering  objects  such  as  non-overlapping 
spheres  or  cylinders  of  infinite  extent  in  the  z-direction  with  circular  or  rectangular  cross  section, 
embedded  in  a  host  medium.  In  these  cases  there  is  a  complete  set  of  functions  to  expand  the 
electric  and  magnetic  currents.  No  restriction  is  imposted  on  the  lattice  two-  or  three- 
dimensional  configuration. 

Both  the  homogeneous  (eigenvalue)  and  inhomogeneous  (plane  wave  source)  problem  were 
formulated.  Two  approaches  were  presented  for  the  three-dimensional  eigenvalue  problem, 
either  by  the  application  of  the  Bloch  theorem  (cf.  Eqs.  (113),  (114))  or  by  solving  the 
homogeneous  algebraic  system  in  Eqs.  (63a),  (63b)  where  the  right  hand  side  is  set  equal  to  zero. 

The  integral  equation  approach  provides  an  easy  transition  from  the  real  to  the  spectral 
domain  either  in  the  two-  or  three-dimentional  lattice  of  dielectric  spheres.  In  Appendix  VII,  the 
proof  that  the  identity  in  Eq.  (99)  is  valid,  is  rather  simple.  But  it  was  only  after  the  integral 

equation  approach  was  applied  that  it  became  evident  that  Z|™  as  given  by  Eqs.  (92),  (1 15)  in 

the  real  domain  is  identical  to  Zj^  as  given  by  Eqs.  (97),  (1 16)  in  the  spectral  domain. 


In  a  lossless  host  medium,  the  sums  in  should  converge  much  faster  in  the  spectral 
domain  rather,  than  in  the  real  domain.  Therefore,  the  computation  of  the  scattered  fields  and  the 
reflection  and  transmission  coefficients  should  be  faster.  In  addition,  as  it  was  shown  in  the 
paper  by  N.  Stefanou,  et.  al.  [6]  (cf.  Table  1  in  their  paper),  a  small  value  of  lmax  (where 
l=l,2,...,lmax)  is  sufficient  for  an  accurate  computation  of  the  scattered  fields.  Consequently,  the 
size  of  the  matrices  involved  in  the  computation  is  small.  Therefore,  there  is  a  small  requirement 
for  memory  as  well  as  computer  time.  Furthermore,  due  to  the  small  size  of  the  matrices,  the 
poor  convergence  problems  for  hard  spheres  in  the  plane-wave  method  or  the  unstable  solutions 
in  the  transfer  matrix  method,  should  not  arise  here.  In  conclusion,  the  approach  presented  here, 
although  limited,  could  be  used  as  a  benchmark  to  compare  and  check  the  results  by  more 
advanced  and  elaborate  approaches  such  as  the  plane  wave  expansion  or  the  transfer  matrix 
approach. 


Appendix  I 

Most  of  the  following  relations  are  taken  from  reference  [19].  The  vector  spherical 
harmonics  are  defined  as  follows: 


Yin(Qr)=  r  J.,/  (rxV)Ylm(^rX 

(I  la) 

• 

YffrtjCU  =  7  l!^-17V(r-|-1Y,m(Qr))> 

[(i+l)(2I+l)]/2 

(Hb) 

Y”"M<nr)=ta^V(r'Y|m<nr))’ 

(lie) 

• 

where  Qr=(0,<p),  1=1 ,2,3, — ,  m—1,-1+1,...,  1-1,1,  and  Ylm(Qr)  are  the 
harmonics.  For  computational  purposes,  the  vector  spherical  harmonics  can 
terms  of  the  3j-symbol  [19]  as  follows: 

scalar  spherical 
be  expressed  in 

• 

Yf,(nr)=  £(2i+iy2(-i)l_m( 1  '  ' 

n=-i  -m  +  p  -pj 

(I2a) 

where  1'  =  1  - 1,1,1  + 1 ,  and 

e±l  =+7j(ex±ey)’ 

(I2b) 

• 

e0  =  ez 

(I2c) 

Here,  ex,  ey,  ez  are  the  unit  vectors  along  the  three  axes  in  cartesian  coordinates. 

The  vector  spherical  harmonics  form  an  orthonormal  set,  i.e., 

• 

Jv,7?r  (Qr)  •  Y^(Qr)dQr  =  5,r5ir8mW. 

(13) 

The  integration  is  over  the  47i-solid  angle.  Moreover,  they  satisfy  the  identities 

er-Y,™(Qr)  =  0, 

(I4a) 

• 

e,-Y,^i.,(nr)  =  -(3jTr)'I'rlra(Qr), 

(I4b) 

«,-Yu-u(£2r)  =  (;rH),2Y,m(Or), 

and 

(I4c) 

• 

er  x  Y,S(nr)  =  j[(^ YM.M(£Jr)+ Y,7_,,(Or )  , 

(I5a) 

er x  Yu+1,1  (^r )  =  jfjiTr)' 2  Ym  (^r )' 

(I5b) 

• 

e,xY,r,,l(nt)=j(^)'2Ylll';(nr). 

(1 5  c) 

With  the  help  of  Eqs.  (Ila)-(Ilc)  and  (I5a)-(I5c)  it  can  be  shown  that: 
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V  X  Y„f  (nr )  =  Xatrf2  (-  Tfe.,1  (Or )  ■ +  ¥*"1-1.1  (Or  )> 


VxYi:i+M(Qr)=j(IiLTf2j±2Y,7;(Qr), 


v  x  ,  (Qr )  =  -  j  iQj#  ¥  Yin  (n  r )  (I« 

The  modified  spherical  Bessel  functions  of  the  first  and  third  kind  [20]  are  defined  as  follows: 


fn/yi 

fi(z)=  V  Wz)> 

V  J  2 


§i(z)=  Tz  Ki+*X(Z)’ 
V2V  /2 


where  z  can  be  a  complex  number.  Notice  that  gi(z)  here  defers  by  a  factor  ^  to  that  defined  in 
Abramowitz.  They  satisfy  the  identities: 

fe~z)fi(z)  =  fi+i(z)’  <I8a) 

(£  +  ^l)f.(z)  =  fM(z),  (I8b) 

fe"l)g|(z)  =  -g|+i(z)’  (I8c) 

(di+ir)gi(z)  =  _g|-i(z)-  (I8d) 

From  these  identities  and  Eqs.  (I6a)-(I6c),  the  following  relations  can  be  shown  to  be  valid: 

Vxf,(OT)YS(Qr)  =  j4(iT)^fl*l(OT')Yu.u(nr)+fe^f.-l<^)V,ru(nr)],  (I9a) 


V  x  fw  (orfe (nr)  =  )'*  f,(or)  Y,I  (£lr  t  <[9b> 

V  X  f,.,(or)Y,;_IJ(nr)=  j°(*)%°r)Ym  (Ori  (I9c) 

v xg|(or)Y,” (Qr)=  -j J gw (ar)Y,j., j («3, )+  (i# g,_, (<rr)Y.;.u (£lr )].  (1 10a) 


V  x  gw  (or)YS„  ,  (Qr )  =  -ja^)^  g.forjY.K  (fl, ), 

V  x  g,_,  My, (nr ) = Mi&f2  g,  (or)YS  (a, ). 

The  outward  radial  Green's  function  in  the  spherical  representation  is  equal  to  [20] 

G(r-r')»  «t?Te"|r"r1  Ef,(<ar<)g|(OT>)Ylln(Dr)Y|,m(nr,) 

1  1  1=0  m=-l 

where  r<  is  the  smallest  of  (r,r')  and  r>  is  the  largest  of  (r,r'). 


(1 10b) 


(1 10c) 
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Appendix  II 


Here,  we  shall  prove  the  relation 

I,„  *  |jm(asinu)epcosuP,m(cosu)smudu  =  2(-l)"'f,(,/if: 


a2  P,m 


/  \ 

p 

Va/p  2-«2y 


(HD 


where  a,p  are  complex  numbers.  The  modified  spherical  Bessel  function  of  the  first  kind  f,(z) 
for  complex  argument  has  been  defined  in  Appendix  I.  The  associated  Legendre  function 

P,m(w)  for  complex  argument  w,  as  well  as  fo  real  argument  x,  is  defined  in  Abramowitz  [21]. 
When  m=l,  we  have  the  relations 

PI‘(cosu)  =  (-l)l(21-l)!!(sinu)1, 


P 


Vp2-<*: 


(  '' 

=  (21  - 1)!!!  -■* 


(II2a) 

(II2b) 


where  (21-1)!!=1.3.5... (21-1). 

We  also  have  the  identity  [22]: 

.  | (sinu)l+l cosh(pcosu)j,(asinu)du  = 
0 

By  setting  u'  =  n  -  u ,  we  obtain  the  relation 


i 

r  \ 

a  f 

1 


a ~ 


(H3) 


J(sinu)l+lepC0SUJ,(asinu)du  =  J(sinu')l+Ie  pC0SU  J|(asinu')du' 

% 

and,  therefore, 

I,,  =  2(— l)1  (21  —  1)!!  |(sinu)l+1cosh(Pcosu)J1(asinu)du  =  2(-l)'f, 

o 

The  relation  (III)  holds  for  m=l. 

When  m=l-l,  we  have  the  relations 

P,,_l(cosu)  =  (-l)M(21-l)!!cosu(sinu)M, 

f  \  f  ^ 

s\-\ 


(H4) 


P 


•  (115) 


■.  1 

-  h\  i)"  p 

f  \ 

a 

Up2-«2  j 

1  Np2-2 

Wp2-«2  J 

(II6a) 

(II6b) 


By  setting  1  ->  1  - 1  in  Eq.  (113)  and  taking  the  partial  derivative  with  respect  to  P,  we  obtain  the 
relation 


(H7) 


cos  u(sin  u)1  sinh(P  cos  u)  J,_,  (a  sin  u)du 


/  \f  \ 

1-1 1 

/  . \ 

(  V 

d 

f  1  /R2  nr2  II  « 

-  f  (i/r2  n2 

a 

ap 

f'-'Wp  “  j 

fi^p  a  V-2 

U2-«2J 

where  we  used  the  identity: 

a^2-^f,-l(z)  =  fl<z).  (118) 

The  change  of  variable  u'  =  n  -  u  leads  to  the  relation 


71 

Jcosu(sin  u)1  efiC0SUJl_1  (a  sin  u)du 

Vi 

Vi 

=  -  |cosu'(sinu')'e  pcosu  JM(asinu')du' 

O 

and,  therefore, 


(H9) 


I,  =  2(— l)1  '(21-1)!!  Jcosu(sinu)1  sinh(Pcosu)J1_,(asinu)du 


=  2(-l)'-1f, 


(IHO) 


We  have  shown  that  the  relation  (III)  holds  for  m=l,  1-1.  Next,  we  shall  show  that,  if  it  is 
true  for  (1-1, m)  and  (1-2, m),  then  it  is  true  to  (l,m).  Therefore,  by  induction,  it  is  true  for  m=l-2, 1- 
1,  ...,  1,0.  For  this  purpose,  we  use  the  identity  [23]: 


P1m(z)  =  ^zP1m1(z)-^P1n:2(z). 


(1111) 


Substitution  of  this  identity  into  Eq.  (Ill)  gives  the  relation 


T  —  21-1  d  t  1+m-l  t 

1l,m  1-m  dp l~m  1-2, m 


-  2(-l)r 


21- 


1-m  dp 


pm 

rl-l 


a 


l+m-l  Dm 
1-m 


1-2 


(w) 


(1112) 


Here,  we  have  assumed  that  Eq.  (Ill)  is  valid  for  (1-1, m)  and  (l-2,m).  After  taking  the  partial 
derivative  and  using  the  identity  (118),  as  well  as,  the  identity: 


(z2  -l)^  =  (l-m  +  l)P,':i (z)-(l  +  l)zP,m(z),  (1113) 


we  obtain  the  relation: 


t  _  2(-l)  If  f-)|  i\  P _ pm 


C-O  +  m-ltC 


(a-O-J—f,.,  (21-l)-nL_Pl“l-(l  +  m-l)Pl^-(l-m)Pl”  , 


where  the  arguments  of  f,,P,m  are  given  above.  The  identity 

(21  -  l)zP,™(z)  =  (1  -  m)P,m(z)  +  (1  +  m  -  l)P,m2(z), 
applied  to  the  relation  above,  leads  to  the  validity  of  Eq.  (III). 

When  we  set  a  =  kpr,(3  =  -y2r,  where  y2  =  [kp  +  a2]^,  Eq.  (Ill)  becomes  equal  to 

71 

Jjm(kpr sin u)e  Y 2 rcos upm (cos u) sjn u(ju 

O 

=  2<-l)'P,ra(j|-)fi(<T2r). 

where 


pm/VL  fcpf  dmP|(w) 
1  Vcti  /  \ct2  /  j,„m 


1  \c2I  V2  /  dwm  w=Y2_’ 

a2 

when  m>0,  and 

r>m (r2 

1  Ict2  /  (|+  mj)  I  Va2  r 

when  m<0.  In  Eqs.  (1115)  we  used  the  identity 

P1m(-w)  =  (-l)'"mP1m(w). 

Finally,  we  use  the  identity 


e-jkp  p  =  |;(-jrj„(kpp)e-jmi''-'''‘». 


(1114) 


(II  14a) 


(1115) 


(II  16a) 


(II  16b) 


(1117) 


(1118) 


where  (p,(p)  and  (kp,  cpu)  are  the  cylindrical  coordinates  of  the  vectors  p  and  kp,  respectively. 
With  the  help  of  this  identity,  of  Eq.  (II 1 5)  and  the  completeness  and  orthonormality  of  the  scalar 
spherical  harmonics  it  is  a  simple  matter  to  prove  the  identity 

e-iVP-v2*  =4ltjr  £  (-l)lf,(a2r)Y,;(nr)Yta(^-,<pl)  (1119) 

1=0  m=-l 

where  a2  is  a  complex  number.  The  complex  scalar  spherical  harmonics  are  defined  as  follows: 

.  (1120) 


(1120) 


II-3 


where  P,m(^-)  is  given  by  Eqs.  (16a),  (16b).  Notice  that  P,m  (z)  is  zero  when  |m|  >  1 .  If  the  real 

part  of  g2  tends  to  zero,  i.e.,  medium  2  is  non-absorptive,  then  Ylm  reduces  to  the  usual  scalar 
spherical  harmonics  and  Eq.  (1119)  reduces  to  the  usual  expansion  of  exp(-jk-r)  in  terms  of 
spherical  harmonics. 


We  define: 


Appendix  III 


f,m(P.N>  =  (-l)m]jm(Xp)e-'|z|Prfe)^ 


where 

Also  we  define 

y  =  [a,2  +  ct2 

1/ 

/2 

where 

?m(P.H)  =  g|(or)P|"(^) 

We  shall  prove  that 

r  =  [p2+z2 

Vi 

flm(P^|z|)=fim(P5N)- 

Here,  a  is  a  complex  number  with  a  positive  real  part,  i.e.,  a^O. 

A)  From  the  identity  (w  a  real  or  complex  number)  [23]: 

C(w)  =  ^[(21  +  l)wP,”(w)-(l  +  m)P,”,(w)} 
it  follows  immediately  that 

W  =  d^rtzi  +  >H#i»  -(1  + 

We  shall  prove  that  fjm  satisfies  this  identity  too.  We  have 


+  i§i(OT) 


1- 


We  use  the  identity  (|x|<l,  real) 

[l  -  x2Jf(x)  =  P,T,(x) +i!±S5raPiTi(x) . 

Then 

-  (21  + 1)±^  fta  - (1  +  =  -(21  +  l)MgJ  P,m  -  (1  +  m)[g,_,  +  ^g, ] 

+  (l-m  +  l)lJFg1P“1. 


(ini) 

(III2) 

(IH3) 

(IH4) 

(III5) 

(III6) 

(III7) 

(III8) 

(III9) 

(III  10) 


We  use  the  identities 


gi-i  +^g,  =-g',  , 


(III  11  a) 


^rgl=gl  +g|+i- 


(III  lib) 


Then 

(21  +  !)(-  -  (I  +  ”)?-,.» 

=  (I  -  m  +  +  g;  [-  (21  +  l)Mp,»  +  (1  -  rn  +  1)C  +  (1  +  m)P,™, } 

From  the  identity  (III6)  it  follows  immediately  that  fjm  satisfies  Eq.  (III7)  as  well. 
B)  From  the  identity  (|x|<l,  real) 

P,m+2  (x)  +  2(m  + 1)-^==  P,m  (x)  +  (1  -  m)(l  +  m  +  l)P,m  (x)  =  0 , 

Vl-x2 

it  follows  immediately  that 


f|,m+2  =  “2(m  +  0“f|,m+l  “  0 ' “  mX1  +  m  +  l)^m  ■ 
We  shall  prove  that  fim  satisfies  this  identity  too. 

From  the  definition  of  fim  we  have: 


(III  12) 


(III  13) 


(III  14) 


(-ir'Nfim*,  =  -]J»*'(Xp4e',|z|hm+1fe)f 

0  (III  15) 


• 

= ]k»,  (Mr* 1  fe)+ iu,  (^)p,ra'fe|-,w  t . 

o 

We  use  the  identities  [23-24]: 

+ 

Then 

4+1  (^p)  =  “(m  +  l)Jm+,  (Xp)  +  XpJm  (Xp), 

(w2  -  l)>lm+,’(w)  -  (m  +  l)wP,m+l(w)  =  Vw2  - 1  P,m+2(w). 

(III  16a) 

(III  16b) 

• 

-  2(m  +  l)|f, =  2(m  +  l)(-l)m  j[jm(Xp)J P,"”1 fe)+i  J™»lP,m*lfe)]e_,W^‘- 

0 

(III  17) 

We  use  the  identities  [23-24]: 

W*p)  =  -UM+t^UiP-p). 

(III  18a) 

• 

Then: 

P,m+2(w)  (1  m)(l  +  m  +  l)P|m  ( w)  -  2(m  +  l)  "  P,m+I(w)- 

V  w  -1 

(III  18b) 

• 

0-  m)0  +  m  +  l)Jm(Xp)P|”fe) 

=  2(ra  +  l)[jm(MfP,mtl(j)+iJm^P)Pr2y  ’ 

(III  19) 

III-2 


Substitution  of  this  identity  into  Eq.  (Ill  17)  leads  to  the  relation  we  want  to  prove,  namely,  f|m 
satisfies  also  the  identity  (III  14). 

C)  We  have  proved  that  both  f,m  and  fjm  satisfy  the  identities  given  by  Eqs.  (III7),  (III  14).  If 
foo,  fn  are  known,  all  the  other  f|m  for  1=1,2, 3,...  and  m=-l,  -1+1,  ...,  1-1, 1  can  be  computed  from 
these  two  identities.  Therefore,  if  we  show  that 

foo(P>M)  =  foo(P>|z|)>  (III20a) 

fn(p,|z|)  =  f\,(p,|z|),  (III20b) 

then  we  have  proven  the  validity  of  Eq.  (III5).  Eq.  (III20a)  is  tabulated  in  I.S.  Gradshteyn  and 
I.M.  Ryzhik  [25].  Notice  that  this  relation  is  valid  only  when  the  real  part  of  a  is  greater  than 
zero.  Eq.  (III20b)  follows  from  Eq.  (III20a)  by  differentiating  the  latter  with  respect  to  p. 

From  the  identity 

e‘jXp=  f]H)mJm(Xp)ejm(<P_<^),  (IH21) 

m~-cc 

we  obtain  the  relation 

2je-jxp+jm<pxd(px  =  27t(-j)m  Jm(^p)ejm(,, .  (III22) 

O 

Combining  Eqs.  (III5),  (III22)  with  the  definitions  of  Y,m(Qr)  and  the  complex  scalar  spherical 
harmonics  Ylm(y/a,cpx)  (cf.  Eq.  1119),  we  obtain  the  following  useful  identity: 


where 


gi(ar)Ylm(Qr)  =  ^  J  ]y,  ^signz)^,^^^1 

—00  — oo 


(III23) 


1  if  z  >  0 


signz  = 


- 1  if  z  <  0 


(III24) 


Appendix  IV 

Here,  we  shall  prove  the  scalar  addition  theorem,  namely,  the  identity; 


g,(a|r-r0|)Ylm(Q  )  =  £  ^Glm;|.m.(r0)fr(ar)Yrm.(Qr) 

"  =0m’=r 


where  a  is  a  complex  number  with  a  positive  real  part,  ro>r,  and 

Gim;rm'(ro)  =  4<-l)l"m2Bim(1'm';rm'')g|(crr0)Y1-j_m-(Qro)  ■ 


(IV 1) 


(IV2) 


I'm' 


The  quantity  B|m(rm';rm")  in  terms  of  the  Wigner  3j-symbol  is  equal  to  [26]: 


Blm(l'm';rm'')  = 


(21+l)(2r+l)(2P+l) 

K.firn 

f  i  r  if> 

4n 

o 

O 

O 

l^m-m"m'; 

(IV3) 


Since  B|m(l'm';rm")  is  different  from  zero  only  when  1'  =  |l-r|,|l-l"|  +  l,...,l  +  l" ,  and 
m'  =  m"  -  m,  the  sum  over  (l',m')  in  Eq.  (IV2)  extends  only  over  these  values. 

When  m>0,  from  the  identity  (III24),  we  obtain  in  operator  form  the  relation  (cf.  Eq.  II  16a) 

g,(°r)Y,m(Or)  =  •  (IV4a) 


where 


P(m)(w)  =  i^iM 

1  dwm 


(IV4b) 


If  we  use  the  usual  addition  theorem  when  ro>r,  namely. 


§o(4'_rol)=47tX  Z  (- I)m" f|”(C5T)gr (OTo )V|*m' (f2r )Y,.  _m. (q,.o  ) ,  (IV5) 


r=0  m’=- 

then  Eq.  (IV4)  yields  the  relation  in  Eq.  (IV 1),  where 


If  we  use  again  the  identity  (III24)  and  perform  the  differentiations  as  prescribed  in  Eq.  (IV6), 
we  obtain  the  relation 


J-m+r-m 


Gim;i'm(«o)  =  47r(-l)l+m'(signz0), 


(IV7) 


The  product  of  the  complex  spherical  harmonics  in  Eq.  (IV7)  can  be  written  as  the  sum  of 
complex  spherical  harmonics,  namely, 


(IV8) 


When  the  real  part  of  a  tends  to  zero,  then  the  complex  spherical  harmonics  become  the  usual 
real  spherical  harmonics  and  B,m(l'm';l"m")  is  given  by  Eq.  (IV3).  The  sum  is  extended  over 
l'  =  |l-r|,|l-r|  +  l,...,l  +  r,  and  m'  =  m"-m.  Only  for  those  values  B|m(l'm';r'm")  is  different 

from  zero.  When  w  =  y/o  is  a  complex  variable,  then  Eq.  (IV8)  expresses  really  an  identity 
between  the  product  of  two  polynomials  and  the  sum  of  a  set  of  polynomials.  The  coefficients 

(-l)m  B,m(rm';lTn")  in  the  sum  are  the  same  for  any  value  of  w.  Therefore,  if  Eq.  (IV8)  is  valid 
for  real  w  less  than  one,  it  is  also  valid  for  complex  w. 

Substitution  of  Eq.  (IV8)  into  Eq.  (IV7)  and  use  of  the  identity  (III24)  yields  the  relation 
5lm;,v('-o)  =  4n(-l)1-”2(sign2o)''"4'''”,*l'‘m'B.m(l'm';rm')g,.(OTo)Y1._in.(nro}  (IV9) 


where  we  used  the  fact  that  m'  =  m"  -  m.  Notice  that  B|m(l'm';l"m")  is  different  from  zero  only 
when  r  =  l-l"  +  2n,  for  1>1",  where  n  =  0,l,...,l",  and  only  when  r  =  l"-l  +  2n,  for  1" >  1 , 
where  n=0,l,...,l.  Therefore,  1  +  1' +  1'  and  m  +  m'  +  nT  are  even  integers  and  Eq.  (IV9) 
becomes  identical  to  Eq.  (IV2). 

The  function  G,m.|W(r0)  can  be  written  in  integral  form  rather  easily.  Multiply  Eq.  (IV 1)  by 

f,*(cn-)Y,tm.(Qr)  and  integrate  over  the  volume  V0  of  the  sphere.  Due  to  the  orthonormality  of 
the  spherical  harmonics,  we  obtain  the  relation 

GtaiW(r0)  =  ^  Jg,(a|r  -  r0|)Ylm(Qr_>0  frMY.Unrjd’x  (IV10) 


where  is  given  by  Eq.  (56). 


IV-2 


Appendix  V 

Here,  we  shall  prove  the  vector  addition  theorem,  namely,  the  identity: 
gi(<+~ro|)V2l(nr-ro) 

=  £  £[0,m:,w(ro)frMY«',(nr) 


(■bWffWYjwfe))] 


1-  a  yy  ZmU'm'VQJy  Vl’V-’1  J  L 

where  a  is  a  complex  number  with  positive  real  part,  r0>r,  and 

^lm;l'm'(**o)  —  7  ,  1/ 

[l(l+l)l(!  +l)]/2 

*  [2a  jnp{J|,+1G  i>m+i;1’,m'+i  (*o )  +  2ajl’-1p[11G|  m-i;1',m'-i(ro)  +  mm'G,m;,w(r0)], 


(V2a) 


V|m;l'm'  (*0  ) =  7 - 21'+'  |7  (“ 1 ) 

[i(i+i>r(i'+i)]/2 


rjf 


*[V2a|"Br_|im-+1(l,-l;rm')Gl  m+l.,'_|  m'+l(r0)- V2P|"B,_|  m-_|(l,l;l'm')Gl  m_,(r0)  (V2b) 

+  m'(l,0;l'm')G|  m;c_i  m(r0)], 

Wln,;,w(r„)=0  (V2c) 

Also,  a["  and  p|"  are  given  by  Eqs.  (37a),  (37b)  while  G,m;iw(r0)  is  given  by  Eq.  (IV2). 

A)  First,  we  shall  compute  U|m;,.m<(r0).  Applying  the  curl  on  Eq.  (VI),  we  obtain  the  relation 
iVx(g,(a|r-r„|)Y,K(nr.,0)) 

=  l[v,m;,vk)fr(^)YlT;(nr)+Ulma.m.(r0)1LVx(f|.(OT)Yff;(Qr)|  (V” 


Use  of  the  identity 


L-F(r)  =  -jr-VxF(r) 


where  L  =  -  jr  x  V  and  F(r)  is  any  vector  function,  together  with  Eq.  (I4a)  leads  to  the  relation 

L  •  (g,  (cr|r  -  rolJYjS  (cir-r0 )) 

=  Z  Ota;nn.(r0X-j)r  •  V  x  (fr  (or  W,  (Or )).  ( V5) 

I'm' 

Using  again  the  identity  (V4)  and  Eq.  (I la)  it  follows  that 

-jr  v(f|.(<n-)Yff;(nr))=-—!-n-L-L(f|.(or)Yl.m.(Qr))=[r(r  +  l)F2f|.(or)Y, v(£lr).  (V6) 

Then,  substitution  of  Eq.  (V6)  into  Eq.  (V5)  yields  the  relation 


V-l 


M&(4-rol)Ym(«r-J) 

=  £  [l'O'  + #  U,m.,w  (r„)fr(<Tr)Y,w  (Qr ) 

I'm' 

On  the  other  hand,  from  the  definition  of  Y,™(Qr)  (of  Eq.  (I la)  and  the  identifies: 

L,Yta(nr)=2al"Y,m„(nr),  (V8a) 

L_Ylm(i2r )  =  2P!"  Y|m_,  (Or  t  (V8b) 

L2Y,m(nr)=rnYlm(Dr),  (V8c) 

where  L±  =  Lx  ±  jLy  and  ajkpj11  are  given  by  Eqs.  (37a)  and  (3 7b), we  obtain  the  relation 

[i(i+i)^Y,;(nr.ro) 

=  kY,,mt,(fi1-ro)+P!nYLm.l(Qr.jK  (V9) 

-  jfct^Yi  m+|(£2r_r0 )—  Pi  Y|  m_[(or_r0)b  +  niY|m(Qr_rf) 

Multiply  Eq.  (V9)  by  g,(cj|r-r0|)  and  apply  the  scalar  addition  theorem,  i.e.,  Eq.  (IV1).  We 
obtain  the  relation 


[l(l  +  l)]^g,  (cr|r  -  r0|)Y™  (nr_ro ) 

=  Xf(X|11^'l.m+l;l,m'(ro)+PrGi,m_i;|'m'(ro)]ex 

I'm' 

-j[arGl.m+i;iw(i'o)-PrGl,m-i;iw(ro)K+mGlm;iw(io)ez)fr(CTr)Y1v(^r) 

Apply  the  operator  L  on  Eq.  (V10)  and  use  Eqs.  (V8a)  -  (V8c).  We  obtain  the  relation 


(V10) 


[l(l  +  l)l^L.(g,(0|r-r„|te(nr.1,)) 

=  yi2"rP"'  +l^l.m+l.l'.m'+l(ro) 

I'm' 

+  2aJ7  _1p[nGl  nV_,  (r0 )  +  mm'G (r0 )]fr  (ar )Y1W  (Qr ) 


(VI 1) 


Eq.  (V2a)  follows  directly  from  Eqs.  (V7),  (VI 1)  and  the  completeness  of  the  spherical 
harmonics. 

B)  Next  we  shall  compute  Vlm;1.m.(r0) .  From  Eqs.  (VI)  and  (V6)  we  obtain  the  relation 

g|(‘+-ro|V'YlnfcVr0) 

IWL  J 

On  the  other  hand,  from  Eq.  (V10)  we  obtain  the  relation 


V-2 


[l(l  +  l)^g,(CT|r-r0|>.YS(nr_ro) 

rm' 

+  P!"G|.m-i;iw(roXx  +  jy)fi'(OT)Y|W(Qr)  +  mG,m;lw(r0)zfr(ar)Y1.m,(Qr)]. 


(VI 3) 


Substitute  the  relations 


x±jy  =  +^frYul(C2r), 
z  =  J¥rYl0(nr) 


(VI  4a) 
(VI  4b) 


into  Eq.  (VI 3)  and  replace  the  products  of  the  spherical  harmonics  by  means  of  Eq.  (IV8)  for 
real  y/o.  Then  Eq.  (VI 3)  becomes  equal  to 


where 


[1(1 + i)l^gi(<y|r  -  r0|V  •  Y,;(nr_ro ) 

-S[a  lmil'm'  (r0  )rf|Vi  (cn-)  +  B  lm;l'm'(r0  >fr_,(ar)]Yrm.(Qr), 

I'm' 

Aim;iv(ro  )=(-0mV¥ 

*  [V2aj"G|  m+|.|'+|  m'+1(r0)Br+1  m-+1(l,-l;l  m  ) 

—  V^P^G,  m_1;|-+,  m-_|(r0)B,.+1  m-_|(l,l;l  -  m  ) 

+  mG|m;i'+i,m'B|'+i,m'(UO;l  m )] 

Blm;1w(r0)  =  (-l)m'V¥ 

*  [V^aj^Gi  m+|;i'_! ,m'+l(ro)B|'-l .m'+l(h— 1;1  m  ) 

-  V2p]nG|  m_1; (*o)B|-_1  m-_,(l,l;l  -  m  ) 

+  mG|m;1'_i  m'Br_!  1  m  )]. 


(VI 5) 


(VI  6a) 


(VI  6b) 


Substitute  the  relations 


arfr+1  (err)  =  arf,'  (err)  -  l'fr  (ar), 
arfr_,(crr)  =  errf,' (ar)+  (l'  +  l)f,.(crr). 


(VI  7a) 
(VI  7b) 


into  Eq.  (VI 5).  Then  from  Eqs.  (VI 2),  (VI 5)  and  the  completeness  of  the  spherical  harmonics, 
we  obtain  the  relation 


(r0)fr(ar)+  B,m;rm,(rn)arf|,(ar)  =  0, 


(VI 8) 


where 


Aw»-(<b)=  Da  + 1)1(1'  +  # Vl„n>o)+  l’Alm;l.„.(r0)-(l'  +  (V19a) 


V-3 


B,m;l-m-(>o)=  [1(1  +  l)f2  Wlm;rm.(r0)-  A,m.l.m.(r0)-Blnl;l.m,(r0).  (VI  9b) 

Equ.  (VI 8)  is  valid  for  all  values  of  r.  Since  there  is  at  least  one  pair  of  values  (zi,  z2),  where 
zi=ari,  Z2=ar2,  for  which  the  determinant 


fl'(Z|)Z|f|'(z,) 

fi’(z2)z2fr(z2l 


is  different  from  zero,  we  conclude  that: 

Alm:l'm'(*b)  —  ®lm;l'm'(**o)  —  ^  ’ 
and,  therefore,  we  have  the  relations 

Vlm;,w(*b)  =  7 - J— n-[-l'A,m;iW(r0)+  (1'  +  l)Blm;lw(r0)], 

W|m;rm'(*"o)  —  ~  ,1/  [Alm;l'm'(lb)~*"  ®lm;l'm'(*b)]" 

[l(l  +  l)]/2 


(V20) 

(V21a) 

(V21b) 


C)  From  Eq.  (I4a)  and  the  identity  V  •  (LY,m)  =  0,  where  L  =  -jr  x  V,  we  obtain  that  relation 

V.(h(ar)Y™(Qr))=0,  (V22) 

for  any  function  h  (ar).  Eq.  (V22)  is  independent  of  the  origin  of  the  coordinate  system,  i.e.,  it  is 
valid  also  for  r-r0.  Therefore,  taking  the  gradient  of  Eq.  (VI),  we  conclude  that 

Y, W,m;,w(r0)V2(f1.(ar)Yrm.(Or))  =  0  (V23) 

1  9 


or 

Ew,m;,w(r0)fi'MYrm.(Qr)  =  0 .  (V24) 

!'m' 

From  the  completeness  of  the  spherical  harmonics  it  follows  that  Eq.  (V2c)is  valid.  Also,  from 
Eq.  (V21b)  we  obtain  the  identity 

Alm:l’m'(ro)+  Blm;lw(ro)  =  0  »  (v25) 

and  Eq.  (V21a)  becomes  equal  to 

U|m;i'm'(ro)~:  21 +1  |y(-l)  B|m.|-m.(r0),  (V26) 

[l(l  +  l)l'(l’+l)]/2 


which  is  identical  to  Eq.  (V2b). 


Appendix  VI 

First,  we  shall  compute  the  Fourier  transform  of  f1(ar)Ylm(Qr)  inside  the  volume  V0  of 
radius  R,  and  the  Fourier  transform  of  g,(ar)Ylm(Qr)  outside  the  volume  V0.  For  this  purpose, 
we  use  the  identity 

eJXt  =  4it £  i;f,0^Km(n,)Ylm(QlL)  (VII) 

1=0  m=-l 

where  k  is  the  magnitude  of  k.  Due  to  the  orthonormality  of  the  spherical  harmonics,  we  obtain 
the  following  two  relations: 

j  f,(OT)Yjnr>JXrd3x  =  4ttYlm(Q  J  Jf,  (cn-lf.O^rVMr,  (VI2a) 

Vo  « 

J  g,(OT)Ylin(Qr)e'x'd3x  =  4JtYln,(£)l)]g,(<3r)f,0>-r>2dr.  (VI2b) 

v«,-Vo  R 

The  two  integrals  over  r  above  are  equal  to 

[f,  (orjf,  (jkr)r2dr  =  Fi  (*-> CT)’  (VI3a) 

J  X  +ct 

0 

F, (X,a)=  crRf.^laR^Om)- jXRf^aR^^OXR)  (VI3b) 

and 

Tgi  (ar)f|(jXr)r2dr  =  -y^-G^a),  (VI4a) 

j  X  +CT 

R 

where 

G, (X,a)  =  aRg1_,(oR)f1(jXR)+jXRg1(aR)fl_1(j?iR)  (VI4b) 

Then,  the  final  result  is  as  follows: 


Jf.larjYjnr^'d’x  =  ^F.fMYjnX),  (VI5a) 

\  g1(<rr)Yln,(Qr)eJ1'd3x  =  ^G,(X,<x)Ylm(nl),  (VI5b) 

\r~  Z  i 


y  =  [K+(y2Y2’ 

(VI6a) 

K  =  +x\  Y2- 

(VI6b) 

where 


Next,  we  shall  compute  the  Fourier  transforms  of  the  vector  functions  in  Eqs.  (55a),  (55b), 
which  we  define  as  follows: 

E£fcta>(X)  =  |E^n)(r)ejVrd3x,  (VI7a) 

v0 

ESmU,)W=  |ESr,(r)eJ’-'dJx,  (VI7b) 

Vcx>-Vo 


where  a=M  or  E.  The  vector  function  E(2^Jn)(r)  is  given  by  Eq.  (10a),  and  E^a,’°ut)(r)  is  given 
by  Eqs.  (40a),  (40b).  Then,  with  the  help  of  Eqs.  (I2a),  (1 10a)  and  the  Fourier  transforms  just 
computed,  it  is  rather  easy  to  show  that 


£(M,in) 

^2,lm 


47iR 


F|(^2)Y|S(£\), 


E 


(M,out) 

2,1m 


4rcR 

^z+Y2 


o,(^.<T2)Y,s(nx), 


(VI8a) 

(VI8b) 


*©">(>.)=■ -f^\  (ikt'2 


t  m 
L  1,1+1 


V21+1 


2 G|_,  (x, CT2  J  ) 


,  (VI8c) 


where 


12  +G22f2 . 


y2  =  k 

Equ.  (IV8c)  can  also  be  written  as  follows: 

Egrw= E^TM*  ^S,(o2R)Vx  x  (f,(jXR)Y;(£lx)). 


(VI9) 

(VI 10) 


Appendix  VII 


We  shall  prove  Eq.  (99),  where  we  set  a2=o,  here.  From  the  integral  form  of  Glm;iw(apq), 

i.e.,  Eq.  (IV10),  and  the  Fourier  transforms  of  g|(ar)Ylm(Qx),f,(cjr)Y|m(Qx),  i.e.,  Eqs.  (VI5a), 
(VI5b),  we  obtain  the  relation 


G 


fjV(M  ClMy* 

hl+v*2  A.2,+v2  ''m' 


(«>.  )Y,m  )e d3X , 


(VIII) 


the  integration  extending  over  the  whole  volume  in  the  spectral  domain.  Substitution  of  this 
relation  into  Eq.  (92)  and  making  use  of  the  transition  relation  given  by  Eq.  (68)  leads  to  the 
validity  of  Eq.  (99).  Notice  that  the  prime  in  Eq.  (92)  can  be  removed  when  the  integral  form  of 

G,m;i'm’(aPq) is  used  (cf.  statement  before  Eq.  (94a)). 
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